"UNCLASSIFIED* 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Pubic  nwortng  burton  for  this  coflectfon  of  information  is  MUnwftBd  to  moqi  1  hour  per  response*  MudnQ  ths  tiros  for  rwriMing  instructions*  starching  aafsflng  dsts  sourest,  gstoehng  and 
maireerrog  thsdato  needed,  and  consisting  and  reviewing  frecolacion  of  rtorroalcn.  Send  oommsrta  mgulng  Msburtan  ssimtos  or  sny  other  taped  of  thto  coiedton  of  information, 
todudnq  suggestions  for  isdudnQ  He  burton*  to  Washington  Hsartfriertari  Services.  Directorate  tor  Intofineaon  Operations  and  Reports*  1215  JtftoraonOevto  Highway*  SUto  1204,  Artnqton. 

VA  222C2-OQ2.  and  to  fts  Oftics  of  Management  andTaudget,  Paperwork  Reduction  Project  (0704-0188).  Washington.  DC  20601 

1 .  AGENCY  USE  ONLY  (Leave  Blank) 

2.  REPORT  DATE 
!April98 

3.  REPORT  TYPE  AND  DATES  COVERED 

Final  !May97  -  !Jan98 

4.  TITLE  AND  SUBTITLE 

Computational  Tools  for  Distributed  Parameter  Control 
Systems 

5.  FUNDING  NUMBERS 

F49620-97-C-0020 

AFRL-SR-BL-TR-98- 

6.  AimOR(S) 

Dr.  Christopher  K.  Allen 

Dr.  Gilmer  Blankenship 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADORESS(ES) 
Techno-Sciences,  Inc. 

10001  Derekwood  Lane,  Suite  204 
Lanham,  MD  20706 


9.  SPCNSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Offfcfe  of  Scientific  Researchy/ffM 
110  Duncan  Avenue,  Rm  B115 
Boiling  AFB 

Washington,  DC  20332-8080 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES 


12a.  DISTRIBUTION/AVAILABIUTY  STATEMENT 

£)i'$'br>  lodtte'i  (Ml 


19980414  070 


13.  ABSTRACT  (Maximum  200  words) 

The  first  phase  of  the  design  and  development  of  a  computer-aided  design  (CAD)  tool  for  distributed 
parameter  control  systems  has  been  finished.  The  kernel  of  the  design  tool  has  been  designed  and 
implemented  for  the  two-dimensional  situation.  It  is  capable  of  modeling  and  simulating  arbitrary  systems 
described  by  partial  differential  equations.  The  implementation  is  based  on  a  finite  element  method  using 
fictitious  domains  and  is  written  in  C++.  As  such,  the  kernel  can  employ  a  regular  grid  for  all  problem 
geometries  and  simulate  moving  and/or  deforming  bodies  without  regridding. 

We  pay  particular  attention  to  the  application  of  viscous,  fluid-flow  control.  The  design  of  an 
incompressible,  viscous,  fluid-flow  simulator  using  the  CAD  kernel  equations  is  presented  in  detail. 
Analytical  work  on  the  control  flow  problem  is  also  presented.  The  work  is  based  upon  perturbation  and 
asymptotic  methods  as  well  as  analysis  using  differential  geometry  and  Mathematica.  The  objective  here  is 
ro  determine  analytic  principles  on  which  to  base  a  feedback  controller  for  flow  control.  The  simulator 
would  then  be  used  to  test  the  designs. 

_  MIC  QTTAT.TT7  INSPECTED  4 


14.  SUBJECT  TERMS 


Distributed  parameter  control,  feedback  control, 
viscous  fluid  flow,  boundary  layer  control 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

UNCLASSIFIED 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UNCLASSIFIED 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

UNCLASSIFIED 


IS.  NUMBER  OF  PAGES 


16.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 


UL 


MSN  75*001  -280-5500 


REF  D 


Standard  Form  298  (Rw.  2-89) 

ftaaafead  by  ANSI  SWL  239-1 S 
299-102 

"UNCLASSIFIED" 


Techno-Sciences,  Inc. 


10001  Derekwood  Lane 
Suite  204 

Lanham,  Maryland  20706 
(301)577-6000 


Final  Technical  Report 


Contract  F49620-C-0020 


-  % 

Computational  Tools  for  Distributed  Parameter 
Control  Systems 


March  31, 1998 

Dr.  Christopher  K.  Allen 
Principal  Investigator 


Abstract 


The  first  phase  of  the  design  and  development  of  a  computer-aided  design  (CAD)  tool 
for  distributed  parameter  control  systems  has  been  finished.  The  kernel  of  the  design  tool 
has  been  designed  and  implemented  for  the  two-dimensional  situation.  It  is  capable  of 
modeling  and  simulating  arbitrary  systems  described  by  partial  differential  equations.  The 
implementation  is  based  on  a  finite  element  method  using  fictitious  domains  and  is  written 
in  C++.  As  such,  the  kernel  can  employ  a  regular  grid  for  all  problem  geometries  and 
simulate  moving  and/or  deforming  bodies  without  regridding. 

We  pay  particular  attention  to  the  application  of  viscous,  fluid-flow  control.  The  design  of 
an  incompressible,  viscous,  fluid-flow  simulator  using  the  CAD  kernel  equations  is 
presents  in  detail.  Analytical  work  on  the  control  flow  problem  is  also  presented.  The 
work  is  based  upon  perturbation  and  asymptotic  methods  as  well  as  analysis  using 
differential  geometry  and  Mathematica.  The  objective  here  is  to  determine  analytic 
principles  on  which  to  base  a  feedback  controller  for  flow  control.  The  simulator  would 
then  be  used  to  test  the  designs. 
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Final  Technical  Report 

Computational  Tools  for  Distributed  Parameter  Control  Systems 


1  Overview 


Under  the  auspices  of  the  U.S.  government's  Small  Business  Innovative  Research  (SBIR) 
program.  Techno-Sciences,  Inc.  has  contracted  with  the  Air  Force  Office  of  Scientific 
Research  (AFOSR)  to  design,  develop  and  build  computational  software  tools  for  the 
simulation  and  design  of  distributed  parameter  control  systems  (Contract  F49620-C- 
0020).  Particular  emphasis  was  to  be  placed  on  applications  of  controlled  fluid  flow.  The 
overall  goal  of  the  project  is  to  develop  a  general-purpose  Computer-Aided  Design  (CAD) 
tcoi  capable  of  modeling  and  simulating  a  wide  range  of  applications  involving  the  control 
cf  distributed  parameter  systems.  This  document  is  a  summary  of  the  work  performed  in 
the  Phase  I  effort. 

We  envisage  the  final  product  of  this  effort  as  a  full-featured,  stand-alone  system  with 
graphical  front  end  similar  to  that  of  existing  modeling  software  systems  such  as 
AutoCAD.  Ultimately,  the  user  would  be  able  to  model  and  simulate  distributed  parameter 
s;.stems  and  design  and  simulate  control  subsystems  for  such  systems.  The  CAD  tool 
would  ailow  the  user  to  easily  build  complex  systems  from  simpler  subsystems  by 
pmviding  him  or  her  with  a  palette  of  such  subsystems  from  which  to  choose.  Thus,  the 
CAD  tool  represents  the  physical  system  as  a  modular  set  of  interacting  subcomponents, 
as  well  as  providing  simulation  and  control  design  capabilities. 

in  Phase  I  we  used  viscous  fluid  flows  for  our  prototype  of  a  distributed  parameter  system. 
We  are  currently  considering  a  two-dimensional,  incompressible,  Newtonian  fluid  in 
amitrary  domains  over  arbitrary  solid  objects.  The  Navier-Stokes  equations  are  used  as 
the  modeling  equations  in  this  situation.  These  equations  are  nonlinear  and  demonstrate 
craotic  behavior  (turbulence).  Moreover,  they  model  certain  physical  behavior  accurately 
enough  to  be  of  great  practical  interest.  Thus,  fluid  flow  control  contains  most  of  the 
features  with  which  we  are  concerned  when  attempting  to  control  a  distributed  parameter 
system. 


We  concentrated  our  efforts  on  the  development  of  a  general-purpose  modeling 
environment  for  distributed  parameter  systems.  In  our  paradigm,  the  objective  system  is 
comprised  of  a  domain  and  a  set  of  boundaries,  some  of  which  may  be  the  surface  of 
ccjects  within  the  domain.  We  also  assume  that  the  system  can  be  modeled  by  a  set  of 
partial  differential  equations  with  boundary  conditions.  The  control  actuators  and  sensors 
are  assumed  to  be  located  on  the  boundaries  in  the  domain.  A  set  of  C++  class  libraries 
has  been  developed  to  model  such  situations.  The  class  libraries  implement  a  finite 
eement/fictitious  domain  representation  of  a  distributed  parameter  system.  Once  the 
target  system  is  defined  using  the  class  library,  it  generates  the  discrete  representation  of 
the  system  suitable  for  numeric  simulation  on  a  computer.  This  representation  is  in  a 
general  format  suitable  for  control  system  analysis.  The  control  actuators  may  control  the 
cependent  variables  or  the  geometry  of  the  boundaries,  while  the  sensors  detect  the  state 


of  the  dependent  variables  at  the  boundary  surfaces.  Initial  tests  of  the  modeling  system 
were  made  for  the  particular  case  of  two-dimensional  Newtonian  fluid  flow. 

Implementation  of  a  finite  element/fictitious  domain  discrete  modeling  system  using  object- 
oriented  programming  was  a  primary  goal  for  us.  This  way  at  the  very  minimum  the 
project  would  yield  a  very  useful  software  library  for  the  general-purpose  modeling  of 
mechanical  and  electromagnetic  systems.  We  initially  attempted  to  build  the  kernel  by 
programming  Matlab,  however,  this  was  found  to  be  unsuitable  for  modeling  the  general 
situation  (the  programming  was  awkward  and  the  problem  setup  times  were  unusually 
long).  The  advantage  of  this  current  software  library  over  other  finite  element 
implementation  is  the  use  of  the  fictitious  domain  formulation  and,  also,  the  object-oriented 
framework.  The  object-oriented  implementation  makes  the  library  modular  and  easy  to 
use,  the  library  being  based  on  an  intuitive  model  of  the  finite  element.  In  fact,  object- 
oriented  programming  is  inherently  designed  to  promote  “reusability”.  The  fictitious 
domain  method  is  a  special  technique  for  handling  the  boundaries  of  the  problem  domain 
and  their  associated  boundary  conditions.  This  method  is  particular  useful  for  treating 
complicated  geometries  and/or  situations  where  the  boundaries  change  with  time.  Thus, 
since  w<=  use  a  fictitious  domain  method  for  simulation,  we  have  the  option  of  actively 
modifying  system  geometries  in  order  to  achieve  control. 

Our  initial  control  objectives  were  to  actively  control  fluid  flow  around  solid  objects  in  order 
to  suppress,  or  at  least  delay,  the  onset  of  turbulence.  Intimately  related  to  this  goal  is  the 
reduction  of  drag  forces  on  the  object.  We  began  studies  on  control  of  both  the 
continuous  Navier-Stokes  equations  and  on  the  discrete  (computer)  approximations  of 
these  equations.  As  of  yet  we  do  not  have  any  conclusive  results  concerning  this  topic. 
However,  our  preliminary  findings  do  suggest  general  principles  for  the  control  of  viscous 
fluids.  In  particular,  we  have  performed  a  boundary  layer  analysis  using  perturbation 
methods  in  order  to  develop  some  physical  principles  on  which  to  base  our  control 
strategy.  The  most  notable  findings  there  include  the  necessity  of  using  multiple  time- 
scales  to  describe  transient  behavior  in  boundary  layers.  Our  current  results  for  both  the 
continuous  and  discrete  control  problems  are  included  in  the  final  section  of  this  report. 

The  work  done  in  Phase  I  may  be  divided  into  three  categories: 

•  Design  of  CAD  Kernel 

•  Viscous  Fluid-Flow  Simulator 

•  Control  Problem  Structure  for  Fluid  Flow 

We  cover  each  of  these  topics  in  a  separate  section  of  this  report. 


2  Design  of  the  CAD  Tool  Kernel 


2.1  Introduction 

The  kernel  of  the  CAD  tool  is  designed  to  model  arbitrary  distributed  parameter  systems 
which  may  be  described  with  a  set  of  Partial  Differential  Equations  (PDEs).  An  initial 
attempt  was  made  to  build  a  prototype  kernel  using  Matlab  and  a  finite  difference 
formulation.  However,  programming  in  Matlab  was  found  to  be  much  too  cumbersome 
and  awkward  when  attempting  to  model  arbitrary  systems.  The  current  prototype  kernel  is 
written  in  C++,  being  based  upon  a  class  library  we  have  written  called  FemLib.  FemLib 


2 


implements  a  Finite 
Element  Method  (FEM) 
for  computer  modeling 
of  distributed  parameter 
systems;  it  is  useful  for 
any  application  using  a 
finite  element 

representation.  The 
implementation  of 
Fern  Lib  is  based  on  a 
high  degree  of 
polymorphism  over  a 
small  set  of  base 
classes.  This  allows 
the  user  to  define  an  •  Figure  1 :  FEM  and  FDM  meshes 

extremely  broad  range 

of  systems  and  system  topologies  using  a  simple  modeling  paradigm.  Once  the  objective 
system  is  described  using  finite  elements,  the  class  library  will  generate  all  the  finite 
element  matrices  and  tensors  (including  fictitious  domain  matrices  -  see  below).  These 
mathematical  objects  provide  the  discrete  representation  of  the  system  which  may  then  be 
manipulated  using  standard  mathematical  techniques.  One  may  employ  Matlab  or  M++  to 
do  these  computations.  Small  class  libraries  are  built  on  top  of  FemLib  to  implement  the 
soecifics  of  the  actual  system  under  study.  To  simplify  matters  and  reduce  development 
time  the  Phase  I  implementation  has  been  done  for  the  two-dimensional  (2D)  situation. 
Therefore,  all  classes  have  been  suffixed  with  the  identifier  "2D"  to  obviate  this  limitation. 
The  three  dimensional  implementation  would  be  directly  analogous. 

Tne  FemLib  class  library  is  capable  of  modeling  solid  objects  embedded  in  the  problem 
domain  using  the  Fictitious  Domain  Method  (FDM),  or  Domain  Embedding  Method  (DEM). 
This  technique  is  typically  used  in  conjunction  with  the  finite  element  method.  The  general 
idea  behind  the  fictitious  domains  is  to  model  solid  bodies  in  the  problem  domain  (usually 
a  connected  subset  of  Ft3  where  d=  2,3)  mathematically  rather  than  using  special  meshes. 
By  a  special  mesh  we  mean  one  that  is  required  to  contour  around  the  object  and,  thus, 
the  object  is  represented  by  the  mesh  itself.  Accordingly,  if  the  object  at  hand  was  to 
change  or  move,  the  mesh  would  have  to  be  regenerated.  In  the  fictitious  domain  method 
a  regular  mesh  modeling  the  problem  domain  is  ordinarily  used.  In  addition  to  the  domain 
mesh,  boundary  meshes  are  superimposed  which  represent  the  surfaces  of  embedded 
objects  in  the  domain.  The  two  types  of  meshes  need  not  conform  at  all,  see  Figure  1  for 
the  example  of  an  infinite  cylinder.  In  the  mathematical  formulation,  dependent  variables 
are  introduced  on  the  superimposed  boundaries  and  are  interpreted  as  Lagrange 
multipliers.  These  multipliers  enforce  the  boundary  conditions  of  the  original  problem. 

The  fictitious  domain  method  makes  the  kernel  implementation  much  easier  and  allows  us 
to  treat  much  more  general  problems.  This  is  because  the  mesh  used  to  model  the 
domain  and  the  mesh  used  to  model  the  embedded  objects  are  independent.  This 
condition  provides  us  with  many  important  capabilities.  It  allows  for  the  modular 
implementation  of  the  CAD  subsystem  palette,  system  components  described  by  a 
boundary  mesh  may  be  added  or  deleted  easily.  Also,  we  are  able  to  model  situations 
which  would  be  very  difficult  to  model,  if  not  impossible,  with  the  finite  element  method 
alone.  For  example,  we  may  treat  rigid  body  motion  within  a  fluid  without  any  regridding  of 
domain.  That  is,  the  boundary  mesh  representing  the  rigid  body  may  move  through  the 
domain  mesh  according  to  the  equations  of  motion.  Body  deformations  within  a  fluid  may 
be  treated  just  as  easily.  For  example,  compliant  membranes  and  flexible  structures  may 
be  modeled  and  simulated  without  expensive  regridding.  Also,  we  may  perform  shape 
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optimizations  on  the  embedded  objects  to  permit  the  maximization  of  certain  design  goals. 
Finally,  we  note  that  the  FDM  allows  for  the  use  of  special  fast  solvers  that  exploit  the 
regular  structure  of  the  domain  mesh. 

The  major  disadvantage  of  the  fictitious  domain  method  is  the  increased  size  in  the 
discretized  system.  Since  the  boundaries  are  represented  as  additional  meshes  we  must 
add  the  new  node  variables  to  the  original  system.  Also,  the  fictitious  domain  method 
requires  an  additional  coding  effort  not  needed  when  using  the  FEM  alone. 

2.2  Finite  Element  Method 

We  implemented  a  standard  finite-element  approximation  scheme  for  partial  differential 
equations,  such  as  that  found  in  references  [1],  [2],  and  [3].  In  this  section  the  particulars 
of  our  implementation  are  covered  to  assist  in  the  use  of  the  FemLib  C++  library.  Due  to 
the  limited  resources  available  for  Phase  I,  the  current  version  of  FemLib  models  the  two- 
dimensional  situation  only.  Therefore,  the  discussion  is  limited  to  this  case.  For  a  detailed 
description  of  finite  element  methods  see  the  above  references. 

-  i 

2.2.1  The  Elements 

The  basic  idea  behind  the  finite  element  method  is  to  decompose  the  problem  domain  Q 
into  a  number  of  smaller  component  regions,  or  finite  elements.  The  solution  of  the  partial 
differential  system  is  then  taken  as  a  piecewise  approximation  across  each  of  the  finite 
elements.  Within  the  element  the  field  variable,  say  u(x),  is  expressed  as  a  linear 
combination  of  interpolation  function  known  as  shape  functions  F,{x),  or  basis  functions, 
associated  with  the  element.  The  set  of  all  shape  functions  for  all  the  elements  of  the 
domain  forms  a  basis  for  a  finite-dimensional  vector  space,  which  we  denote  Vh(d)  (h  is  a 
parameter  we  include  to  emphasis  the  geometric  size  of  the  individual  elements).  It  is  this 
function  space  Vh(Q)  where  we  take  the  approximate  solution  to  the  partial  differential 
system.  The 

approximate  solution 
will  satisfy  the  integral, 
or  weak,  form  of  the 
partial  differential 
equations  where, 
typically,  we  use  the 
shape  functions 
themselves  as  the 
testing  functions 
(Galerkin’s  method). 

Part  of  the  art  of  finite 
elements  is  choosing 
the  element  shapes 
and  the  shape 
functions  so  that  the 
function  space  l/^Q)  •  Figure  2:  Finite  element  in  grid  and  with  local 

contains  a  good  coordinates 

approximation  to  the 
original  problem. 

Although  FemLib  is  setup  to  use  any  geometric  shape,  such  as  triangles  or  octagons,  we 
shall  concentrate  solely  on  quadrilateral  elements.  This  simplifies  the  discussion  since  the 
shape  functions  are  easily  described  in  the  local  coordinate  system.  Because  we  are 
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using  the  fictitious  domain  method,  regular  grids  are  typically  used  and  this  restriction 
presents  no  shortcoming  to  the  modeling  process.  An  example  situation  is  shown  in 
Figure  2.  In  the  figure  a  trapezoidal  region  in  Cartesian  (x,y)  space  is  decomposed  into 
four  quadrilateral  elements.  Also  depicted  in  the  figure  is  a  single  element  shown  in  its 
local  coordinate  system.  For  the  quadrilateral  elements  this  is  a  Cartesian  system  for 
which  we  give  the  coordinates  (r,s).  Both  r  and  s  take  values  in  the  closed  interval  [-1,1]; 
that  is,  in  the  local  coordinate  system  all  the  elements  are  squares!  To  convert  back  to  the 
original  problem  domain,  there  exists  for  each  element  a  mapping  <p (r,s)  from  the  local  (r,s) 
coordinates  to  the  global  ( x,y)  coordinates.  This  mapping  is  how  the  geometry  of  the 
system  is  captured  in  the  discrete  model.  We  shall  see  that  the  mapping  is  completely 
defined  in  terms  of  the  vertices  and  the  shape  functions  of  the  finite  elements. 


2.2.2  The  Nodes 


An  important  part  of  the  finite  element  method  is  the  use  of  node  locations  in  the  domain 
q.  These  are  points  in  Q  usually  having  special  significance  with  respect  to  the  elements. 
For  example,  they  tend  to  be  vertices,  centroids,  and  edge  bisectors.  This  situation  is 
depicted’in  Figure  2.  There  are  nine  nodes  in  the  figure  which  are  the  vertices  of  the 
elements.  Each  node  has  an  index  /  which  is  its  index  into  the  discrete  system  (in  the 
figure  the  indices  run  from  0  to  8).  This  index  does  not  change  throughout  the  simulation. 
(Each  node  may  have  multiple  element  indices,  however.  See  the  discussion  concerning 
Figure  3.)  Thus,  it  is  possible  to  change  the  structure  of  the  finite  element  matrices  simply 
by  relabeling  the  nodes’  system  indices.  Note  that  the  nodes  may  be  shared  by  more  that 
one  element,  this  is  how  the  topology  of  the  system  is  captured  in  the  discrete  system;  it 
results  in  coupling  between  the  nodes. 


c) 


d) 


The  number  and 
position  of  the  nodes 
on  an  element 
depends  upon  the 
interpolation  scheme 
used  in  the  simulation. 

Typically,  a  polynomial 
approximation  is  used 
and,  therefore,  only  the 
order  of  the 

approximation  is 
variable.  Figure  3 
shows  the  node 

locations  for  the  four 
interpolation  schemes 
implemented  in 
Fern  Lib.  Note  that  as 

the  order  of 

approximation 

increases  the  number  of  nodes  in  the  element  must  increase.  Also  shown  in  the  figure  are 
the  element  indices  n  for  the  nodes.  Each  finite  element  maintains  a  dynamic  array  of 
nodes  that  it  contains  and  the  element  index  is  the  means  by  which  it  references  its  nodes. 
This  index  is  different  from  the  more  ubiquitous  system  index  in  that  each  node  may  have 
a  different  element  index  depending  upon  which  element  is  referencing  it.  The  node 
locations,  as  indicated  in  the  figure,  have  local  coordinate  values  of  either  -1 , 0,  or  +1.  For 
example,  nodes  0, 8,  and  2  of  the  biquadratic  element  in  Figure  3  d)  have  (r,s)  coordinates 
(-1,-1),  (0,0),  and  (+1,+1),  respectively. 


>  Figure  3:  Node  positions  and  indexes:  a)  Constant,  b) 
Bilinear,  c)  Biquadratic  serendipity,  and  d)  Biquadratic 
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2-2.3  The  Shape  Functions 


Here  we  list  the  shape  function  for  each  different  element.  Defining  rn  and  sn  as 


we  have  the  following: 
For  the  constant  element 


F0(r,s)  = 


if  (r,s)e[— l,+l]x[— 1,+1] 
if  (r,s)  € [-1,+1] x  [-1,+1] 


For  the  bilinear  element 


F„{r,s)=^-r„r^-sns)  for«  =  0,...,3  (3) 


For  the  biquadratic  serendipity 

Fn(r,  s)  =  ^  (l  -  rnr\ 1  -  sns\r„r  +  sns  - 1)  for  n  =  0, ...  ,3 

=  ^-(i— r2X1-5) 

F5(r,5)=i(l  +  r)(l-52)  (4) 

F6(r,5)=^(l-r2)(l  +  5) 

^(^^(l-rXi-s2) 

For  the  full  biquadratic 

F„  ( r ,  s)  =  ^  (l  -  rnr\ 1  -  sns)rnrs„s  for  n  =  0, ...  ,3 
F4(r,5)  =  -|(l-r2)(l-s> 

F5(7-,5)  =  +^(1  +  r)(l-j2)- 

;  (5) 

T’6(r,s)=+-(l-'-2)(l  +  .s> 

F7(r,s)  =  -±{l-rll-s2} 

W,s)  =  -±(l-r2ll-s2) 

Note  that  each  shape  function  Fn  is  unity  at  the  node  n  and  zero  at  all  the  other  nodes  on 
the  element.  That  is,  in  local  coordinates  we  have  the  following  relation: 
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(6) 


j  1  if  m-n 

F“(r"’5j=  0  if  m*n 


where  (rm,sn)  are  the  local  coordinates  of  element  node  m.  In  the  finite  element  method 
the  values  of  field  variable  u  are  computed  only  at  the  discrete  node  locations.  Thus,  the 
shape  functions  provide  the  interpolation  of  u  between  these  node  locations. 


Z2A  Integration 

The  finite  element  method,  as  opposed  to  the  finite  difference  method,  is  based  upon 
integration.  Therefore,  we  need  to  know  how  to  integrate  functions  over  the  elements,  this 
is  where  the  mapping  (p  comes  into  play.  Recall  from  Figure  2  that  (p  maps  the  local 
element  coordinates  (r,s)  to  the  global  coordinates  {x,y)  in  the  domain  where  the  element 
is  located.  It  can  be  varified  that  this  mapping  is  given  by 


(x,y)  =  <p(r,s)  =  Y,(x„,yn)Fn(r,s)  (7) 

-  » 

where  the  summation  is  taken  over  all  the  nodes  n  of  the  element  and  the  pairs  (xn,yn)  are 
the  global  coordinates  of  the  n,h  node.  Say  we  wish  to  compute  the  integral  of  some 
function  f(x.y)  over  some  finite  element  T.  From  the  calculus  we  know  that  this  can  be 
done  in  the  local  coordinates  (r,s)  of  T using  the  change  of  variable  formula 

jj  /('V,  y)  dxdy  =  J  J  f{(p(r,  drds  (8) 

where  <3cp (r,s)/d(r,s)  is  the  Jacobian  of  the  transform  <p(r,s).  In  FemLib,  all  the  finite  element 
integrals  are  computed  using  this  formula.  The  integrals  are  computed  numerically  using 
various  quadrature  schemes  that  may  be  selected  by  the  user  (e.g.,  Simpson,  Bode, 
Gaussian). 

In  practice,  the  shape  functions  are  actually  associated  with  the  nodes  of  the  system  and 
are  indexed  by  the  system  index  i.  Thus,  the  shape  function  F-,  is  defined  piecewise  over 
each  element  that  contains  node  i.  The  shape  function  used  on  each  element  depends 
upon  the  particular  element  index  n  of  the  node,  however,  the  system  is  setup  in  such  a 
way  that  the  shape  functions  are  continuous.  Note  also  that  each  shape  function  F)  has 
compact  support  since  it  is  nonzero  only  on  the  elements  that  contain  node  /.  In  this 
context  the  set  of  shape  functions  {F{x,y)}  for  a  basis  for  vector  space,  the  vector  space  of 
approximation  functions.  Since 


Ft(Xj,yj)  = 


Jl  for. j-i 
[0  for  j  *  i 


(9) 


where  (xj,  yj)e R2  is  the  coordinate  position  of  node  j,  any  approximation  of  the  field 
variable  ufxy)  is  expanded  in  terms  of  these  basis  functions  according  to 

u(x,y)«J^uiFi(x,y)  (10) 


where  the  u ,•  are  (the  expansion  coefficients)  the  values  of  u(x,y)  at  the  node  locations 
(Xi,y,).  It  is  the  goal  of  the  finite  element  method  to  find  the  set  coefficients  {u}. 
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2.3  Fictitious  Domain  Method 


2.3.1 


To  see  the  main  idea  of  the  fictitious  domain  method,  consider  the  example  of  a  classical 
Dirichlet  problem,  as  in  reference  [4].  The  problem  is 


au  +  vV2m  =  /  in© 
u  —  g  ony 


(ID 


where  a  and  v  are 
positive  real  parameters, 
u  is  the  dependent 
variable,  and  f 
represents  external 
forces.  The  boundary 
conditions  are  given  by 
the  function  g  which  is 
prescribed  along  the 
boundary  of  the  domain 
Y=5oo.  This  situation  is 
depicted  graphically  in 
Figure  4  where  the 
original  domain  <n  is 
embedded  in  a  larger 
domain  Q  whose 
boundary  is  T. 

The  Variational  Formulation 


•  Figure  4:  Dirichlet  Problem  Domain  and  Boundary 


In  the  spirit  of  finite  element  method,  we  want  to  express  the  above  problem  in  a 
variational,  or  weak,  form.  This  is  done  by  multiplying  the  original  equation  by  a  testing 
function  v{x)  then  integrating  the  Laplacian  term  by  parts.  Modulo  technical  details,  the 
resulting  equation  should  hold  for  any  such  testing  function  v.  By  properly  choosing  the 
set  of  allowed  testing  functions  we  can  eliminate  the  boundary  term  arising  from  the 
integration  by  parts  to  secure  the  following  description: 


Find  u  in  Vg  such  that 

J  [c£«v  + vVm  •  Vv]rfQ  =  J/vJQ  VveF#  (12) 

n  n 

where  Vg  and  Vo  are  the  function  subspaces  given  by 

Vg=\teH\(o)\u  =  gony},  ^ 

F0=/f‘(©)  =  ^e/T1(©)|«=Oonr}. 


This  is  the  problem  that  the  standard  finite  element  method  approximates.  The  field 
variable  u  is  expanded  in  terms  of  the  shape  functions,  according  to  Eq.  (10),  where  the 
node  values  u,  along  the  boundary  y  are  given  by  g  (either  directly  or  with  some  sort  of 
integral  averaging).  This  ensures  that  the  approximation  for  u  is  in  the  set  Vg.  Therefore, 
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we  need  only  determine  the  value  of  u,  for  each  internal  node  of  the  domain  ©.  Note  that 
the  set  of  shape  functions  associated  with  these  internal  nodes  form  a  basis  for  a  finite¬ 
dimensional  space  contained  in  y0- '  Thus,  these  shape  functions  form  a  set  of 
independent  testing  functions,  one  for  each  internal  node.  Applying  Eq.  (12)  to  each  of 
these  functions  yields  a  linear  system  of  N  equation  in  N  unknowns,  where  N  is  the 
number  of  internal  nodes  in  co.  That  system  may  then  be  solved  using  standard  numeric 
techniques  to  find  the  values  of  u,for  the  internal  nodes  and,  thus,  the  full  approximation. 

2&2  The  Fictitious  Domain  Formulation 

Solving  the  variational  form  using  the  finite  element  method  requires  that  the  domain  © 
and  boundary  y  are  approximated  by  the  finite  element  grid.  If  the  boundary  is 
complicated  this  may  be  a  difficult,  or  expansive  task.  Or,  if  the  boundary  changes  with 
time  the  domain  must  be  regridded  at  each  time  step.  The  fictitious  domain  method 
sidesteps  this  issue  by  embedding  ©  in  a  larger,  simpler,  domain  Q  then  enforcing  the 
boundary  conditions  on  y  mathematically. 

The  fictitious  domain  method  may  be  formalized  using  the  Lagrangian  for  the  original 
system.  This  procedure,  in  effect,  results  in  a  minimum  energy  argument.  The  formal 
Lagrangian  L  for  Eq.  (1 1 )  is  given  by  the  following: 

L(u)  =  —  j|aw2  +  v|Vw|2Jf©  -^fudw  V«  e  Vg  (14) 

^  <o  <0 


It  is  well  known  that  minimizing  the  above  Lagrangian  with  respect  to  u  in  V9  yields  the 
weak  solution  to  the  Dirichlet  problem  (indeed,  the  necessary  condition  dL/du=0  produces 
the  variational  formulation  for  the  problem).  Note,  however,  that  this  minimization  may  be 
interpreted  as  a  constrained  minimization  due  to  the  requirements  of  the  set  Vg, 
specifically,  u  must  equal  g  on  the  boundary  y. 

Using  the  Lagrange  multiplier  theorem  we  may  reformulate  the  weak  problem  as  the 
following  unconstrained  minimization 


mm  max 


Z,(m,A)  =  +  v|Vu|2  dco-^fudco+^X(u-g)dy 


VueV, l  e/Tl/2(y) 


(15) 


where  this  time  the  solution  space  V  may  be  chosen  much  more  arbitrarily.  (We  take 
A  s  H~]  ~(y) ,  the  dual  of  the  trace  space  Hy2(j)  .  simply  out  of  convenience.  In 
actuality  we  assume  that  X  is  in  the  larger  space  L2( y). )  In  fact,  we  shall  choose  V so  as  to 
expand  the  actual  problem  domain  to  that  of  Q.  Taking  the  functional  derivative  of  L(u,X) 
with  respect  to  both  u  and  X  yields  the  variational  formulation  for  the  "fictitious  domain" 
Dirichlet  problem 


Find  u  in  i/and  X  in  HU2(y)  such  that 


|  [auv  +  vVii  •  Vv]rfQ  =  J  JvdQ.  +  J  Xvdy 

n  n  r 

\  ,u(u  -  g)dy  =0 


VveV 

(16) 

V^e/r,/2(y) 
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where  we  take 


V  =  H'0(  Q).  (17) 


After  finding  the  solution  uto  this  problem  on  the  larger  domain  Q,  the  restriction  u|B  will  be 
the  weak  solution  to  the  original  Dirichlet  problem. 

The  boundary  conditions  are  enforced  in  Eq.  (16)  by  the  term  involving  the  Lagrange 
multiplier  X.  We  must  solve  for  the  new  variable  X,  as  well  as  for  the  dependent  variable  u, 
in  order  to  solve  the  entire  system.  Therefore,  mathematically  the  situation  is  somewhat 
more  complicated  but  this  presents  no  real  problem  when  solving  the  system  numerically. 
By  comparing  the  two  terms  on  the  RHS  of  the  first  of  Eqs.  (16)  we  see  that  the  Lagrange 
multiplier  X  appears  as  a  fictitious  external  force,  or  impulsive  forcing  function,  acting  along 
the  boundary  y.  This  force  is  exactly  that  which  would  be  caused  by  the  original  boundary 
condition^  acting  on  u.  For  our  current  system  that  force  is 

X  =  (18) 

dn 

where  n  is  the  unit  normal  to  the  boundary  surface.  This  may  be  shown  using  Green’s 
theorem  applied  to  the  variational  formulation  in  (16). 

Note  that  fictitious  domain  formulation  is  readily  solved  (numerically)  using  a  finite  element 
approach.  Indeed,  one  may  choose  the  auxiliary  domain  Q  as  a  regular  subset  of  Rtf  that 
conforms  to  the  coordinates  (e.g.,  a  coordinate  rectangle  for  cartesian  coordinates).  We 
may  then  divide  Q  into  regular,  homogeneous  finite  elements  so  that  the  resulting  matrices 
will  demonstrate  a  high  degree  of  structure  and  symmetry. 

2A  Computer  Implementation 

Here  we  describe  the  computer  software  architecture  for  implementing  the  CAD  tool.  An 
initial  investigation  using  a  finite  differencing  technique  programmed  in  Matlab  was 
attempted  early  in  the  project.  It  was  found  that  the  performance  of  Matlab  was 
prohibitively  slow  when  doing  the  setup  calculations  (computer  modeling)  and  successive 
over-relaxations  for  the  time-stepping  solutions.  Formulating  the  programming  paradigm 
for  system  modeling  was  also  very  clumsy  in  the  Matlab  programming  language. 
Therefore  we  have  implemented  the  problem  formulation  in  the  C++  programming 
language  using  the  M++  class  library  for  matrix  operations.  From  the  C++  code  we  can 
create  MEX  files  for  use  with  Matlab.  In  this  manner  the  resulting  system  matrices  can  be 
send  to  Matlab  to  do  the  command,  numeric,  and  graphical  processing.  Since  matrix 
operations  and  manipulations  in  Matlab  are  relatively  fast,  and  we  have  a  large  variety  of 
visualization  tools  available  for  viewing  the  results,  this  seems  to  be  a  reasonable  setup. 

We  also  describe  here  some  of  the  finite  element/fictitious  domain  modeling  strategy  since 
it  is  relevant  to  the  implementation  and  to  the  structure  of  the  control  problem.  For 
example,  in  order  to  circumvent  the  well-known  divergence  instability  condition,  which 
causes  wild  pressure  fluctuations  in  the  numeric  integration  of  the  Navier-Stokes 
equations,  we  provide  the  capability  of  using  two  different  grids  for  the  pressure  and 
velocity  variables.  Thus,  the  implementation  allows  for  different  sets  of  basis  function 
expansions  on  the  same  grid.  This  leads  to  different  (discrete)  system  structure. 
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2A1  The  Finite  Element/Fictitious  Domain  Library  (FemLib) 

The  computer  modeling  portion  of  the  project  has  been  implemented  as  a  stand  alone, 
general  purpose  set  of  C++  class  libraries.  The  most  basic  library  is  the  finite  element 
library,  called  FemLib,  that  implements  all  the  requirements  of  a  finite  element/fictitious 
domain  formulation.  This  library  is  useful  for  any  application  using  a  finite  element 
representation  of  a  distributed  parameter  system.  The  implementation  is  based  on  a  high 
degree  of  polymorphism  over  a  small  set  of  basic  classes.  This  allows  the  user  to  define 
an  extremely  broad  range  of  systems  and  system  topologies  using  a  simple  modeling 
paradigm.  Once  the  objective  system  is  described  using  finite  elements,  the  class  library 
will  generate  all  the  finite  element  (including  fictitious  domain)  matrices  and  tensors  with 
the  M++  class  library  representation.  These  mathematical  objects  provide  the  discrete 
representation  of  the  system  which  may  then  be  manipulated  using  standard 
mathematical  techniques.  One  may  employ  Matlab  or  M++  itself  to  do  these 
computations.  The  other  class  library  on  which  we  shall  focus  is  the  Navier-Stokes  class 
library.  This  is  a  small  set  of  classes  that  are  used  when  the  governing  system  equations 
are  the  Navier-Stokes  equations. 

-  ) 

To  simplify  matters  and  reduce  development  time  the  Phase  I  implementation  has  been 
done  for  the  two-dimensional  (2D)  situation.  Therefore,  all  classes  have  been  suffixed 
with  the  identifier  ”2D"  to  obviate  this  condition.  The  three  dimensional  implementation 
would  be  directly  analogous. 

At  the  heart  of  the  finite  element  implementation  are  the  virtual  base  classes  FeElem2D, 
FeComplex2D,  and  FeGrid2D.  Also  important  is  the  stand-alone  class  FeNode2D.  A 
class  hierarchy  diagram  is  shown  in  Figure  5.  The  standard  convention  is  to  represent 
inheritance  via  a  directed  arrow  from  the  derived  class  (subclass)  to  the  base  class 
(superclass).  From  this  diagram  we  see  that  the  class  FeElem2D  is  the  base  for  all  finite 
element  classes.  As  it  stands,  it  is  simply  a  container  of  FeNode2D  objects,  which 
represent  the  nodes  within  the  finite  element.  From  there  the  derived  classes  add  the 
various  attributes  and  features  particular  to  any  finite  element. 

The  FeComplex2D  class  is  a  container  class  for  FeElem2D  objects.  This  makes  sense 
since  a  finite  element  grid  is  simply  a  collection  of  finite  elements  that  fit  together  "nicely1  (a 
complex)  such  that  the  topology  is  easily  inferred.  The  class  FeComplex2D  is  a  base  for 
two  types  of  sub-complexes,  one  defining  the  problem  domain,  FeDomain2D,  and  one 
defining  the  problem  boundaries,  FeBoundary2D.  Also,  FeBoundary2D  has  the  derived 
class  FeEmbObj2D  representing  embedded  objects  in  the  fluid  domain.  Therefore,  in  the 
current  two-dimensional  implementation  FeDomain2D  contains  two-dimensional  finite 
elements  (derived  from  FeBody2D)  and  FeBoundary2D  contains  one-dimensional  finite 
elements  (derived  from  FeFace2D).  Finally,  the  complete  finite  element  grid  is 
represented  by  the  FeGrid2D  base  class.  This  class  has  three  major  attributes.  The  first 
is  a  FeDomain2D  object  representing  the  domain  of  the  grid  and  the  second  is  a  container 
of  FeBoundary2D  objects  which  represent  all  the  boundaries  in  the  grid's  domain.  In  this 
manner  we  may  associate  a  separate  boundary  object  for  each  different  boundary 
condition.  The  third  is  a  container  of  embedded  objects.  A  separate  FeEmbObj2D  object 
is  contained  for  each  embedded  solid  in  the  domain. 


11 


•  Figure  5:  Finite  Element  Class  Library  Inheritance  Diagram 

Notice  that  there  are  two  base  classes  of  FeElem2D,  they  are  FeBody2D  and  FeFace2D. 
The  FeBody2D  class  designates  finite  elements  containing  area  (volume  for  the  3D 
situation)  and  representing  the  actual  component  bodies  from  which  the  grid  is  composed. 
The  FeFace2D  class  is  meant  to  represent  components  of  boundaries  within  the  domain. 
Since  there  are  two  types  of  boundaries  in  the  domain  that  we  wish  to  consider,  element 
body  faces  and  free  boundaries  (i.e.,  boundaries  of  embedded  objects),  there  are  two 
subclasses  derived  from  FeFace2D.  The  first,  FeBodyFace2D,  represents  the  face 
boundary  of  a  FeBody2D  object  while  the  second,  FeFreeFace2D,  represents  an  arbitrary 
boundary  in  the  domain.  This  situation  is  natural  because  the  shape  functions  of  the  body 
faces  are  dependent  upon  the  interpolation  scheme  used  in  the  body  element.  The 
interpolation  scheme  and,  therefore,  the  shape  functions  for  a  free  boundary  are  arbitrary. 
However,  we  are  currently  revising  this  implementation  so  that  the  FeFace2D  objects 
naturally  know  whether  or  not  FeElem2D  ” 

they  are  body  faces.  When  Cardinal  nNodes;  //  no.  of  nodes 

completed  the  new  Cardinal  n Vertices;  // no.  of  vertices 

implementation  will  eliminate  ArrayOfR2s  arrVertices;  // element  vertices 

the  need  for  the  two  derived  ArrayOfFeNodes  arrNodes;  //  element  nodes 
classes. 

FeBody2D 

We  consider  the  FeElem2D  Interpolation  {Const,  BiLinear,  BiQuadratic, ...}; 
hierarchy  in  more  detail.  virtual  Real  ShapeFunc(lndex  iNode,  R2  pt)  =  0; 

Figure  6  shows  code  virtual  R2  ShapeGrad(lndex  iNode,  R2  pt)  =  0; 

excerpts  from  these  finite 
element  classes.  In  the  FeFace2D 

figure  the  most  important  Interpolation  {Const,  Linear,  Quadratic, ...}; 
attributes  and  member  of  virtual  Real  ShapeFunc(lndex  iNode,  Real  pt)  -  0; 

each  of  the  base  classes  are  virtual  R2  ShapeGrad(lndex  iNode,  Real  pt)  =  0, 

listed.  (The  type  Cardinal  L  - 

represents  cardinal  numbers,  1  •  Figure  6:  Attributes  of  the  Finite  Element  Classes 

the  type  Real  represents 
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elements  of  the  set  of  real  numbers  R,  while  the  type  (class)  R2  represents  elements  of 
ff.)  The  FeElem2D  class  consists  primarily  of  an  array  of  vertices  (points  in  F?)  and  an 
array  of  FeNode2D  objects.  The  objects  are  contained  in  arrays  since  the  vertices  and 
nodes  of  a  finite  element  are  indexed  according  to  the  type  of  element.  The  FeBody2D 
and  FeFace2D  classes  define  an  additional  attribute,  the  enumeration  Interpolation.  This 
attribute  determines  the  type  of  interpolation  scheme  used  by  the  element,  specifically,  the 
degree  of  polynomial  representing  the  shape  function.  The  interpolation  scheme  for  each 
element  can  be  changed  by  simply  setting  the  appropriate  flag.  Two  member  function  are 
also  declared  in  each  class,  ShapeFuncQ  and  ShapeGradQ.  These  functions  represent, 
respectively,  the  shape  function  and  the  its  gradient  for  the  finite  element.  The  signatures 
of  the  functions  in  each  of  the  two  classes  are  different  since  the  local  coordinates  of  a 
(2D)  body  are  two-dimensional  while  the  local  coordinates  of  a  (2D)  face  are  one¬ 
dimensional.  The  two  functions  are  declared  as  pure  virtual  functions  and  must  be  defined 
in  any  instantiable  base  class,  since  a  finite  element  without  a  shape  function  is  useless. 
The  functions  are  meant  as  placeholders  so  one  can  access  the  shape  function  and  its 
gradient  without  knowing  the  particulars  of  the  finite  element  itself. 

Finally,  the  classes  FeQuad2D ,  FeRect2D,  FeTrian2D,  FeLine2D  represent  actual, 
instantiable  finite  elements.  The  class  FeQuad2D  represents  an  arbitrary  quadrilateral 
body.  The  FeRect2D  represents  a  rectangular  body  and  is  derived  from  FeQuad2D 
(since  a  rectangle  is  a  quadrilateral).  Due  to  the  regular  structure  of  a  rectangle  some  of 
the  coordinate  functions  and  shape  functions  are  faster  to  compute  than  with  a 
FeQuad2D.  The  FeTrian2D  class  implements  an  arbitrary  triangular  body  while  the 
FeLine2D  represents  a  boundary  element  comprised  of  a  line  segment. 

Now  consider  the  FeGrid2D  hierarchy.  As  we  have  mentioned,  the  base  class  FeGrid2D 
contains  three  major  attributes,  a  FeDomain2D  object,  a  list  of  FeBoundary2D  objects, 
and  a  list  of  FeEmbObj2D  objects.  The  topologies  of  these  contained  objects  are 
determined  by  the  sharing  of  nodes  between  the  finite  elements.  That  is,  the  finite 
elements  recognize  their  neighboring  elements  by  sharing  common  nodes,  typically  at 
vertices.  The  FeGrid2D  class  is  simply  a  container  class  of  finite  element  objects  and  no 
provisions  are  made  to  create  any  such  grids.  That  task  is  accomplished  through  the  use 
of  derived  classes.  There  are  currently  two  derived  classes  of  FeGrid2D,  the  class 
FeGridRect2D  which  implements  a  rectangular  grid  and  FeGridFile2D  that  implements  a 
grid  defined  by  a  formatted  input  file  (currently  unfinished).  Since  we  are  using  rectangular 
grids  in  the  test  bed  problems  it  was  convenient  to  provide  a  special  implementation  for 
that  situation.  For  more  arbitrary  situations  the  user  can  define  a  grid  with  an  ASCII  file 
containing  the  nodes  and  elements,  and  their  positions  and  connections. 

The  last  class  we  mention  in  the  finite  element  class  library  is  Felntegrate2D.  This  class 
computes  the  discrete  system  data  from  FeGrid2D  objects.  That  is,  it  computes  the 
matrices  and  tensors  common  to  a  Galerkin-type  finite  element/fictitious  domain 
approximation  of  the  model  problem.  The  discrete  representations  of  the  bilinear  and 
triiinear  forms  appearing  in  the  variational  formulation  are  computed  in  terms  of  the  basis 
of  shape  functions  over  the  elements.  For  example,  the  standard  inner  product  of 
functions  over  the  grid  space  has  matrix  entries  (the  basis  matrix)  given  by 

Ay=\FiFJdQ  (19) 
n 

where  F,  is  the  shape  function  for  the  (internal)  node  i,  while  the  H'( Q)  inner  product  has 
matrix  entries  (the  stiffness  matrix)  are  given  by 
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Cy  =  J"  VFj  ■  VFj  dQ. .  (20) 


2.4.2  The  Navier-Stokes  Class  Library  (FluidLib) 

The  viscous  fluid-flow  simulator  is  implemented  as  a  small  C++  class  library,  FluidLib, 
which  sits  on  top  of  FemLib.  FluidLib  currently  realizes  the  finite  element/fictitious  domain 
formulation  of  the  2D  incompressible  Navier-Stokes  equations.  It  also  provides  for  the 
simulation  of  boundary  sensors  and  actuators  for  a  feedback  control  system. 


•  Figure  7:  FluidLib  classes 

The  simulator  functions  by  first  retrieving  from  FemLib  the  matrices  and  tensors 
representing  the  linear,  bilinear,  and  trilinear  forms  specific  to  the  integral  formulation  of  the 
Navier-Stokes  equations  (we  use  a  primitive  variable  formulation).  These  matrices  and 
tensors  are  assembled  into  a  nonlinear,  discrete  system  that  approximates  the  time 
evolution  of  a  viscous,  incompressible  fluid.  These  equations  are  integrated  in  time  using 
a  three-step  "0-scheme".  This  method  is  a  generalization  of  the  Peaceman-Rachford 
scheme  and  was  introduced  by  Glowinski  in  the  mid  1980's. 
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Currently,  there  are  four  classes  in  FluidLib,  NsBndForms2D,  NsDomForms2D, 
NsSystem2D,  and  NsTimeStep2D.  The  relationships  between  these  classes  are  shown 
in  Figure  7.  In  the  figure,  an  oval  represents  a  class  where  a  rounded  rectangle 
represents  an  actual  object  of  an  instantiable  class.  The  classes  NsBndForms2D  and 
NsDomForms2D  are  essentially  helper  classes  that  compute  the  bilinear  and  trilinear 
forms  associated  with  the  weak  form  of  the  Navier-Stokes  equations.  The  class 
NsSystem2D  actually  holds  these  forms,  represented  as  matrices  and  rank-three  tensors 
(these  data  are  stored  as  M++  matrices  and  arrays,  respectively).  NsSystem2D  also 
owns  the  pressure  and  velocity  grid  (FeGrid2D)  objects.  The  time  integration  of  the 
Navier-Stokes  equations  is  done  in  the  class  NsTimeStep2D.  This  class  owns  the 
NsSystem2D  object  and  maintains  the  dependent  variables  (i.e.,  the  velocities  and 
pressure).  As  mentioned  above,  this  class  implements  the  0-scheme  for  numerical 
integration.  For  the  mathematical  details  of  the  fluid-flow  simulator  see  the  next  section. 


3  Viscous  Fluid-Flow  Simulator 


Here  we’describe  the  techniques  used  to  build  the  viscous,  fluid-flow  simulator.  Most  of 
the  material  discussed  is  implemented  in  the  C++  library  FluidLib.  The  presentation  is 
rather  technical  since  we  discuss  the  details  encountered  in  the  implementation.  There 
would  be  no  loss  of  continuity  if  the  reader  were  to  proceed  directly  to  Section  4. 

3.1  Mathematical  Formulation 

The  governing  equations  are  the  incompressible  Navier-Stokes  equations.  The  situation 
we  describe  is  uniform  flow  around  a  solid  object  embedded  in  the  fluid  domain.  The 
situation  is  described  mathematically  by  the  following  partial  differential  system: 


U  +  (V  •  U)U  +  Vp 

-  vV2U  =f  inQ 

V  •  U(r)  =  0 

inQ 

U(r)  =  Umf 

onr0 

dn 

onf! 

+"1 

Xj 

II 

onrc 

(21) 


where  U(f)eRd  (d=2,3)  is 
the  velocity  vector,  p(t)<=  R 
is  the  pressure 
(normalized  to  fluid 
density),  veR+  is  the 
viscosity,  f(f)eRd  is  the 
external  force  on  the  fluid, 

Uj,/eRd  is  the  fluid  velocity 
at  infinity  (i.e.,  uniform  flow 
velocity),  c(f)eRd  is  the 
control,  neRd  is  the  normal 
vector  to  the  respective 
boundary,  Q  is  the 
problem  domain,  r0  is  the 
Dirichlet  boundary,  I")  is  •  Figure  8:  Fluid  flow  problem 

the  free-flow  boundary, 
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and  rc  is  the  control  boundary  (boundary  of  the  embedded  solid  object).  See  Figure  8  for 
a  graphical  depiction  of  this  situation.  Note  that  rc  might  be  a  function  of  time. 

The  function  space  of  solutions  for  velocity  U (t)  will  be  denoted  V(Q)  (or  simply  V)  while 
the  solution  space  for  pressure  p(t)  will  be  P(Q)  (or  simply  F).  We  will  also  consider 
(vector-valued)  functions  spaces  on  the  boundaries,  denoted  by  S(r0)  and  6(rc).  Note 
that  V(Q)  should  be  a  subspace  of  the  Sobelov  space  H\q.).  For  convenience  of 
reference,  Table  1  collects  the  definitions  of  these  spaces  and  other  function  spaces  used 
in  this  formulation. 

Now  consider  the  weak,  or  "variational",  formulation  for  the  system.  Using  a  fictitious 
domain  formulation  (see  Section  2.3.2)  for  the  internal  boundary  we  have 


Ju 

Jn 

■VdQ  + 

L<°- 

V)U  •  VdQ 

-LfV-VJC J  +  v| 

VU • VVdQ  = 

J 

[  f  •  V  dd  +  f  A  • 
In  Jrc 

VdT 

VVeF0(Q) 

tv 

Jo 

•  VqdQ  ■ 

=  0 

V9s/>(Q) 

f  u 

•b  dT  = 

f  uinf 

•b  dr 

Vb  e  5(T0) 

Jr0 

lr„ 

j," 

■b^T  = 

f  c-bdT 

lrc 

Vb  s  B(TC) 

(22) 


Note  that  the  boundary  conditions  on  n  are  absorbed  naturally  into  the  weak  formulation 
(this  is  by  design).  Referring  to  the  first  line  in  the  above,  the  "testing"  functions  V  live  in 
V0{Q),  which  is  the  subspace  of  V(Q)  given  by 

F0(Q)  =  {veF(Q)|v|ro=o}  (23) 


That  is,  Vyo)  is  the  subspace  of  functions  in  V(Q)  which  are  zero  along  the  boundary  r0. 
A  final  point  with  respect  to  the  first  line  in  (22)  is  that  the  vector  function  A  is  the  Lagrange 
multiplier  resulting  from  the  fictitious  domain  formulation.  It  enforces  the  internal  boundary 
conditions  on  rc.  However,  this  function  is  now  an  additional  dependent  variable  for  which 
to  solve. 

In  the  second  line  of  (22)  we  take  our  testing  functions  from  the  subspace  P1(Q)cP(Q). 
This  is  the  subspace  of  pressure  functions  that  are  zero  along  the  boundary  rt  (rather 
than  r0)  or,  more  precisely, 

;>(«)  =  {?e/>(0)|?  |r,  =  0 }  (24) 

The  motivation  for  this  particular  choice  of  test  functions  will  become  clear  later.  The 
function  spaces  are  summarized  in  the  following  table: 


Space 

Description 

m 

Solution  space  for  U(Q,  generally  a  subspace  of  Sobolev  space  [/-/(Q)]1^ 

Vo(Q) 

Space  of  velocity  testing  functions,  elements  of  V{Q)  which  are  0  on  r0 

m) 

Solution  space  for  p (/),  generally  a  subspace  of  Sobolev  space  W’( Q) 

Space  of  pressure  testing  functions,  elements  of  P(Q)  which  are  0  on  r. 

SCTo) 

Space  of  boundary  condition  testing  functions,  subspace  of  [Hu'\r0)]a 

B(  rc) 

Solution  and  testing  space  for  A(0,  subspace  of  Sobolev  space  [H'l£(rc)]° 

Table  t :  Function  spaces  for  Navier-Stokes  fictitious  domain  formulation 


3.2  System  Discretization  -  Finite  Elements  and  Fictitious  Domains 

For  the  discrete  approximation  we  use  a  finite  element  approach.  Therefore,  a  finite 
triangulation  3*  of  Q  is  assumed  (the  "triangulation"  does  not  necessarily  consist  of 
triangles,  just  some  collection  of  polytopes).  This  triangulation  depends  upon  the 
parameter  h,  the  triangulation  becoming  increasingly  finer  as  h  is  approaches  zero.  Next 
a  function  space  l/*  is  formed  such  that  each  element  therein  may  be  represented  by  a 
polynomial  of  some  finite  degree  n  over  each  triangle  T  in  the  triangulation  3ft.  This 
condition  may  be  written  more  precisely  as 

Vh (Q)  =  jv,  e  K(Q)|v*|r  e  Poly(n)yT  e  3, }  (25) 

where  Poly(n)  is  the  space  of  polynomials  over  Q  with  degree  less  than  or  equal  to  n.  Also 
associated  with  the  triangulation  is  a  set  of  node  locations  {x}h  on  the  domain  Q.  These 
nodes  are  the  sample  points  for  the  dependent  variables,  that  is,  their  values  are  defined 
at  these  discrete  locations. 

-  f 

Incidentally,  the  bulk  of  the  programming  effort  lies  in  the  triangulation  of  the  domain  Q  and 
the  representation  of  the  function  spaces  over  the  triangulation.  Once  the  triangulation 
and  function  spaces  have  been  defined,  the  resulting  discrete  system  may  be  attacked 
using  existing  software,  such  as  Matlab  or  M++. 

3.2.1  Shape  Functions 

Obviously  V(Q)  has  a  finite  basis.  The  basis  we  shall  consider  is  the  set  of  standard 
shape  functions  (F(x)}  for  finite  elements  in  the  triangulation.  These  shape  functions  have 
the  following  properties 


F^j)  = 


for  j  =  i 
for  j  *  i 


(26) 


where  x;eRd  is  the  coordinate  position  of  node  j.  Moreover,  F(x)  has  support  only  on  the 
finite  elements  Te 3h  to  which  node  /  is  common.  In  other  words,  F(x)  is  nonzero  only  on 
the  elements  7e3h  in  which  point  x,  is  contained.  We  shall  denote  the  set  of  triangles 
containing  xi  as  3/^i)  so  this  property  may  be  formally  written 

F,(x)  =  0  Vx  €  3A  \  3;,(0  (27) 

Thus,  we  see  that  the  shape  functions  provide  an  interpolation  scheme  between  nodes. 
Note  that  we  may  have  different  sets  of  shape  functions  for  the  same  triangulation 
depending  upon  the  degree  n  of  polynomial  space  l/(Q). 

Since  the  shape  functions  form  a  basis  for  ^(Q)  we  may  approximate  the  dependent 
variable  Ue  \/'(Q)  as  a  finite  expansion  of  these  functions.  Henceforth,  for  concreteness 
we  assume  that  the  fluid  domain  Q  is  two  dimensional  and  we  separate  the  vector  U=(u,v) 
into  its  scalar  components  to  save  storage  space  in  the  resulting  discrete  system.  Then 
the  basis  function  expansion  for  U  is 
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(28) 


fX“<-  Fi(xA 


u  = 


Likewise,  we  may  invoke  a  similar  procedure  for  the  other  dependent  variables.  Letting 
{G/x)}  be  the  set  of  shape  functions  for  the  pressure  approximation  space  F*(Q)  and 
{H{x)}  be  the  set  of  shape  functions  for  the  boundary  function  approximation  space  6^(rc), 
we  have 


and 


-  * 


p  =  2>;Gf(x)  (29) 


V  J  ) 


where  A=(X,,p). 

3.2.2  Index  Sets 


Inserting  the  above  expansions  into  the  weak  formulation  (22)  and  using  the  shape 
functions  themselves  as  the  testing  functions  (yielding  Galerkin's  method)  produces  the 
system  of  discrete  matrix-vector  (and  tensor)  equations.  Before  we  do  this,  however, 
some  indexing  sets  are  introduced.  For  the  set  of  sample  nodes  {x}h  on  the  function 
space  l/’fQ)  let 


/  =  { 0,1,2,  ...,Ny  - 1 }  (31) 

be  their  index  set,  where  Nv  is  the  number  of  nodes  in  {x}h.  Now  let  /oc/  be  the  set 

/0.{ieJ|Jft*)|r#=o}  (32) 

which  is  the  set  indices  for  the  "internal  nodes"  not  on  the  Dirichlet  boundary  r0.  Thus, 

Vq  (Q)  =  span  {Ft  (x)  |i  e  J0  }  .  (33) 

(The  space  V0h(Q)  is  the  approximation  of  V0(Q.)).  We  also  define  the  set  of  nodes  indices 
on  the  boundary  r0  as 


IB=I\I0.  (34) 

Regarding  the  nodes  of  the  other  approximation  spaces,  we  define  their  index  sets  in  a 
similar  fashion.  The  index  set  for  (G/x)}  is  given  by 

J  =  { 0, 1, 2 . tfp-l}  (35) 
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where  NP  is  the  number  of  nodes  in  the  space  P^Q).  The  index  subset  for  P^h(Q)  is  given 

by..  ,  .  _  ^ 

7,  =  {fe7|Gi(*)|ri  =0}  (36)  ' 

which  is  the  set  of  “internal  pressure  nodes”  not  on  the  free-flow  boundary  IY  Thus, 

P\  (H)  =  span  fa  (jc)  |i  e  7,  } .  (37) 

The  index  set  for  the  set  of  pressure  nodes  on  the  boundary  r,  may  be  identified  as 

JBmJ\Jx.  (38) 

(Note  the  notational  difference  between  lB  and  JB.  The  elements  of  lB  are  indices  of  nodes 
on  r0  and  the  elements  of  JB  are  indices  of  nodes  on  IV)  Finally,  the  index  set  for  the 
basis  {Hffl}  is  defined 

AT=  {  0,1,2,  (39) 

where  NB  is  the  number  of  nodes  in  the  boundary  space  B^Tc).  The  following  table 
summarizes  the  various  node  index  sets. 


Set 

Decription 

/ 

Indices  of  all  nodes  in  the  velocity  grid 

b 

Indices  of  velocity  nodes  not  on  the  boundary  r0;  “internal  velocity  nodes” 

Ib 

Indices  of  velocity  nodes  on  the  boundary  To 

J 

Indices  of  all  nodes  in  the  pressure  grid 

Ji 

Indices  of  pressure  nodes  not  on  the  boundary  Tu  “internal  pressure  nodes” 

Jb 

Indices  of  pressure  nodes  on  the  boundary 

K 

Indices  of  all  nodes  on  the  boundary  rc 

•  Table  2:  Node  index  set  descriptions 


These  sets  allow  the  immediate  description  of  all  the  approximation  spaces  for  the  field 
variables.  These  function  spaces  are  collected  below  in  Table  3  for  convenience. 
Compare  these  spaces  with  those  of  Table  1 . 


Space 

Description 

l//,(Q)=span  {F{x)\iel}  c  V(Q) 

Approximation  space  for  U (t) 

t/on(Q)=span  {F{x)\lel0}  c  V0(Q) 

Space  of  velocity  testing  functions 

Pft(Q)=span  [G{x)\ieJ\  c  P( Q) 

Approximation  space  for  p(t) 

P1n(Q)= span  (G/x)|/e71}  c  P,(Q) 

Space  of  pressure  testing  functions 

BT(Tc)  =span  {H(x)\ieK\cz  B(TC) 

Approximation  and  testing  space  for  A (t) 

•  Table  3:  Approximation  function  spaces  for  Navier-Stokes  problem 
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3.2.3  Discrete  Equations 


Now  we  insert  the  function  expansions  of  (28),  (29),  and  (30)  into  the  weak  formulation 
(22)  to  yield  a  set  of  ordinary  differential  equations.  Using  (F(x)xO,  OxF/x)|/,/e/o}  as  the  set 
of  velocity  testing  functions  (for  the  first  line  of  Eqs.  (22))  we  get 


[4y]u°+v[5,y]u°  +  ]u)u  +  ([Cp ]v)u  - [Dy  ]p 1  -[Ey]\  =  f1  +  W*(ufl,vfl) 
[4>.]v°+v[5iy]v°  +([C^]u)v  +  ([C|t]v)v-[Aj]p1  -[^]P  =  f'  +  WK(us,  v5) 

[Dy  ]Tu°  +[Dy  ]Tv°  =0 

[£,yfu°=c  (40) 

[Fy]Tv°=d 


where  we  have  employed  the  following  definitions  for  the  discrete  dependent  variables: 


u  =(«,•  \i  e/0  u  IB  )  u°  =(m,- |/ e/0  )  nB  =[ut\i  el B)  P  ( Pj  \j  e  A  )  X-[Xk\keK ) 
v=(v,-|f  e/0u/s)  v°  =(v,|f  e/0)  \B=(vi\ieIB)  ps  ■  [pj\j  ^JB)  ^-{^k\keK) 

(41) 


The  forcing  vectors  are  defined  as 


“inf  =(uinf(*/)-a*  I*  e/a)  fS  =( \jX  F‘  e/°] 

V  inf  =  (  Ujnf  (xi)'ay  I1  e/fl)  f-v=^J"  fy  F.  dQ\i  e/0  j 


Jr  cH,dT\i  eK 
dsf  Jr  dH;dT\ieK 


(42) 


Note  that  we  have  used  a  simple  pointwise  approximation  scheme  for  the  boundary 
values  on  r0.  The  node  values  for  U(x,)=[u(x,),v(xl)]  are  specified  directly  by  the  equations 
uB=uint  and  vB=vw  rather  than  using  an  averaging  (integral)  approach.  (Actually,  this  is 
equivalent  to  using  Dirac  delta  functions  as  the  testing  functions.)  This  makes  the  system 
simpler  and  does  not  introduce  any  additional  error  (this  error  has  the  same  order  as  that 
already  incurred  by  the  shape  function  interpolation  scheme). 

Note  that  the  advection  components  in  (40)  may  be  decomposed  into  the  internal  node 
and  boundary  node  components  as  well.  The  result  is  given  by 

(tq£]»)u  +  ([C^]v)u  =  ([C^]u°)u°  +  ([C$]u°)ufl  +  ([q£]ufl)u° + ([q*]u*)ufl 

+  ([^]V°>°  +([C£]v<y  +([C!k]vBy  +  ([C^]vfl)u* 
([C*]u)v  +  ([C^vjv  =  ([C*]u°)v0  +  ([C*]u°)vfl  +  ([C*]ufl)v0  +  ([C*]u*)v5 

+  ([Cjiyy  +([Cj>°)vs  +  ([C;jfc]vfl)v0  +  ([C^]vfl)vs 

The  matrices  and  tensor  elements  in  (40)  are  given  by  the  values 
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4-  =  j>V 

B>~~1 


da 


M.^1.  &l  dh 

dx  dx  +  dy  dy 


D; 


xt  dlLc. 

1  Jn  dx  J 


D: 


,J  Jn  dx 
dE 


da 


v  f  dr  i 

t  =  -r-Gjda 
v  J  ady  J 


for  ij  e/ 

ci  -J 

r  SFj 

F  ~r~  Fk  da 
a  ‘  dx  k 

for  i,j  s/ 

da  for/,  /  e/ 

•.'•••••'I 

•  8F: 

II 

Fi-^LFkda 
a  1  dy  1 

for  Q  el 

for  /  el0,j  €  J 

for/  el0,j  eJ 

FjHjdT  for  /  el,j  eK 

(43) 


The  functions  Wx  and  WY  in  the  LHS  of  (40)  represent  the  contributions  due  solely  to  the 
boundary  nodes.  The  values  of  these  expressions  are  known  functions  of  uB  and  v. 
They  are  defined  by 


WY(us,vs)  =  -[4y]uS -[By] us  -Rf]PS  -{([c^]uS)uS  +([CiH-]vS)ufl}  where  /  el0,j  eIB 

WK(uVV=-^.]vs -[By]**  -[^]P®  -{([Cp]us)vs  +([C^]vfl)vs}  where/  el0,j  eIB 

(44) 

where  the  terms  in  braces  are  only  added  if  the  advection  terms  are  decomposed  into  the 
boundary  node  and  internal  node  contributions.  The  conditions  at  the  end  of  the 
definitions  are  meant  to  indicate  that  the  matrix  and  tensor  elements  in  the  above 
equations  have  the  same  definitions  as  those  in  (43)  only  that  the  indices  belong  in  the 
index  set  lB  rather  than  /0.  In  essence,  we  are  just  separating  a  larger  system  built  from  / 
into  the  unknown  components  k  and  the  known  components  fe. 


3.3  Numerical  Integration:  Operator-Splitting  Methods 

Now  that  we  have  extracted  a  discrete  approximation  to  the  continuous  system  we  may 
use  established  mathematical  techniques  to  solve  the  system.  We  employ  an  operator¬ 
splitting  technique  for  the  time  integration  that  decomposes  the  original  system  into  linear 
and  nonlinear  subsystems.  In  this  manner  we  may  associate  the  (difficult) 
incompressibility  condition  with  the  linear  system  and  leave  the  nonlinear  system  free  (i.e., 
unconstrained).  A  conjugate-gradient  algorithm  solves  the  resulting  linear  Stokes  problem 
and  a  least-squares  algorithm  may  be  used  to  solve  the  nonlinear  advection  problem. 
Here  the  general  properties  of  using  operator-splitting  methods  for  numerical  integration 
are  covered.  In  the  following  section  we  apply  the  techniques  to  the  discretized  Navier- 
Stokes  system. 


3.3.1  Operator-Splitting  Methods 

Operator  splitting  methods  are  a  family  of  techniques  for  numerically  integrating  general 
operator  equations.  They  are  based  on  the  notion  of  decomposing  an  (possibly  nonlinear) 
operator  into  separate  components  which  are  individually  simpler  than  the  original.  To 
demonstrate  this  numerical  integration  technique,  we  shall  use  the  following  system  as  a 
model.  Consider  the  operator  equation  on  the  Hilbert  space  H 


0  +  A  0  =  0 
^(O)  =  0o 


(45) 
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where  A:H->H  is  a  (possibly  nonlinear)  operator  on  H,  <p(\)eH,  and  the  over  dot  indicates 
time  differentiation.  Now  suppose  A  may  be  nontrivially  decomposed  as 

A  =  A,  +  A2.  (46) 

The  above  system  is  to  be  integrated  using  operator  splitting.  A  straightforward  technique 
is  the  Peaceman-Rachford  scheme  [5]  which  we  shall  describe.  First,  however,  we  must 
introduce  the  discrete  time  and  its  notation. 

A  finite-length  time  step  Af  >0  is  chosen.  The  dependent  variable  <p  is  then  indexed  for  the 
discrete  times  nAf  according  to  the  following: 

*,■#»*).  (47) 


NOTE: 

Another  potential  technique  for  solving  such  problems  is  the  use  of  semigroup  methods. 
This  idea  is  based  on  the  generalized  solution  of  the  operator  equation  0  =  A(j>  which  is 
represented  abstractly  as  <p(t)  =  eA' (f)(0) .  These  techniques  have  been  well  studied  for 
linear  systems  [10], [11]. 

3.3.2  The  Peaceman-Rachford  Scheme 

With  the  Peaceman-Rachford  scheme  we  divide  the  time  interval  [nAf, (/Hi  )A/j  into  two 
subintervals  [nAf,(n+1/2)A/]  and  [(n+1/2)Af,(n+1)A/j.  Assuming  that  <pn  is  known,  we 
approximate  the  value  of  <Pn*v2  using  a  backward  Euler  step  with  respect  to  An  and  a 
forward  Euler  step  with  respect  to  A2.  Then  the  value  of  <p^  is  approximated  from  cp^v2 
using  the  same  process  only  with  the  roles  of  At  and  A2  reversed.  The  result  is  that 


K+Vl  ~ (h 

At/2 


+  A  (<f>n+l/2)  +  A(4n)  ~  0’ 


<h+ 1  ~  ^n+1/2 
At/2 


+  At(<f>n+l/2)  +  A2  (tfrn+l )  -  0, 


(48) 


must  be  solved  successively  for  (pn±v2  and  tp^.  Note  that  the  above  is  an  implicit  scheme. 
However,  the  idea  here  is  to  choose  At  and  A2  so  that  the  nonlinear  and  linear 
components  are  separated  making  the  solution  of  the  above  easier. 

We  examine  the  stability  and  accuracy  of  the  Peaceman-Rachford  scheme  by  considering 
a  simple  example.  Let  H=  Rw,  <peR  ,  and  A=[Aj]  be  a  positive-definite,  symmetric  N*N 
matrix.  The  solution  of  Eq.  (45)  is  then  given  by 


m  =  e-lAlj]l<t>0  (49) 


which  may  be  decomposed  into  modal  equations  using  the  eigenvalues,  {A.,},  and 
eigenvectors,  [ej,  of  [>y.  The  solution  in  terms  of  the  component  modes  is  given  by 

N 

(f>(t)  =  ^ 0' (t)  where  <t>‘(t)  =  '$5  i  =  l,..., N.  (50) 

i=i 
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The  quantity  (p\t)  is  the  modal  component  of  cp(f)  for  the  eigenvalue  X,  (i.e.,  the  projection  of 
<p(tj  onto  e,}  and  <po  is  the  component  of  <po  in  the  e,  direction.  By  linearity  we  are  able  to 
track  the  modes  separately.  To  apply  the  Peaceman-Rachford  method  we  impose  the 
following  decomposition  of  the  matrix  [>4J: 

[AiJ]  =  a[AiJ]  +  p[AiJ]  (51) 

where  a,pe(0,1 )  with  a+p=1 .  Solving  Eqs.  (48)  for  cp^  with  A=a[4)]  and  A2=^[A^  yields 

0n+I  (52> 

(6, -is  the  Kronecker  delta  function)  so  that  the  discrete  approximation  to  Eq.  (49)  is  given 
by 


0n=([s*]fy/3K])  ([5,y]+y«(4-])  (to/l-y/w)  - (53) 

Doing  a  modal  expansion  for  this  discrete  solution  yields  the  following  solution  in  terms  of 
the  modal  components  %  of  %\ 


0i  = 


1  At  1 
1 - OU: 

2  ' 

,  Al  1 
1  -i CL X: 

v  2  ' 


lyA 

iyA 

v  2  j 


b‘n  for  i  =  l,...,iV. 


(54) 


Since  |(1-x)/(1+x)|<1  for  all  x>0,  the  above  equation  shows  that  for  all  discrete 

times  ri> 0.  This  fact  implies  the  stability  of  the  Peaceman-Rachford  method,  at  least  for 
linear  systems.  The  above  equation  also  implies  linw<p„'=0  for  all  /,  which  we  expect 
from  Eq.  (49)  and  the  positive  definiteness  of  [Ail.  To  determine  the  accuracy,  set  ja=A ft, 
and  introduce  the  rational  function  R(p)  defined  by 


*00- 


a 

1- ~rju 
a 

1  +  TM 


V 


1  +  f" 


(55) 


Talyor  expanding  fl(n)  about  the  point  n=0  yields 


1  ,  a2+ap  +  p  ,  4x 

^(/i)  =  l-/i  +  -Ai‘ - - 2 - A+0(p)  (56) 


which  when  compared  to  the  Taylor  expansion  of  e1*  about  n=0, 
*-"=1  (57) 

2  o 


shows  that  we  have  at  least  second-order  accuracy  for  any  values  of  a  and  p.  The  best 
values  of  a  and  p  conforming  to  the  constraint  a+p=1  are  a=p=1/2.  This  gives  us  "almost" 
third  order  accuracy  since  (a2+ap+p2)/4=3/1 6=1/6+1/48. 
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The  definition  of  R(\x)  also  points  out  the  main  disadvantage  of  the  Peaceman-Rachford 
method,  namely,  it  converges  slowly  for  stiff  systems.  For  large  values  of  p,  and 
consequently  for  large  values  of  Atk„  the  value  of  R{ n)  approaches  unity.  This  implies  that 
<p„'  converges  slowly  to  zero  as  o-» oo.  From  Eq.  (50)  we  know  that  the  true  solution  <p\ t) 
will  converge  rapidly  to  zero  as  f— *».  Therefore,  the  Peaceman-Rachford  method  is  not 
able  to  capture  fast  transients  unless  At  is  taken  to  be  extremely  small.  This  condition  is 
also  true  of  the  Crank-Nicolson  integration  method. 

3.3.3  The  0-Scheme 

The  0-scheme  is  a  generalization  of  the  Peaceman-Rachford  scheme  and  was  introduced 
by  Glowinski  in  the  mid-1980s  [12].  In  this  method  the  time  interval  [nAt, (rn-l)Afl  is  broken 
into  three  subintervals.  The  result  is  an  unconditionally  stable  integration  scheme  which  is 
able  to  capture  the  behavior  of  stiff  systems  much  more  accurately  when  using  larger 
time-steps  At. 

Again  we  use  Eq.  (45)  as  our  model  system.  The  0-scheme  begins  by  picking  a 
paramfctlr  06(0,1/2)  and  a  time  step  At  >0.  Assuming  that  <pn  is  known,  we  then  solve  the 
following  system,  consecutively,  for  <jve,  (pmis,  and 


+  4(<U)  +  4(«»i-.>  -»■  (58) 

To  check  the  performance  of  this  algorithm  with  respect  to  the  Peaceman-Rachford 
method  again  consider  the  simple  situation  of  H=  Rw,  (peR'v,  and  A=[Aj|  a  positive-definite, 
symmetric  NxN  matrix.  We  proceed  as  before,  solving  Eqs.  (58)  for  cp^  with  A|=a[4J  and 
4>=p[4]  where  <x,pe(0,1 ),  a+p=1 .  Setting  0=1  -20,  the  solution  is 

<t>n+ ,  =  ([5y  ]  +  a6At[Ay  ])"2  ([5ff  ]  -  /J0A t[Av ])'  ([5,y ]  +  pff  A t[Av ])  '  ([5,y  ]  -  ad'  At[Ai} ]}t>n  (59) 

and  the  corresponding  discrete  modal  decomposition  is 

(t-/^.)2*  =  (60) 

(1  +  aOA ai)2n  (1  +  pd'At^f 

To  explore  the  stability  for  stiff  systems  define  the  rational  function  R(y) 


R'(fi)  = 


(l-^)2  (1  -affn) 
(1  +  aO/j.)2  (1  +  p0'n) 


(61) 


Since 


a 


(62) 
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requiring  a>p  should  insure  better  convergence  than  the  Peaceman-Rachford  scheme 
when  integrating  stiff  systems.  In  fact,  the  larger  the  value  of  a  compared  to  p  the  better 
the  convergence  properties.  However,  we  should  also  consider  the  accuracy  of  the  0- 
scheme.  To  do  this  Taylor  expand  ff(p)  about  n=0  to  obtain 

R'( n)  =  l-n+  j[i  +  (P2  -  a2 )( 262  -46  +  l)]#i2  +  0(fi3) .  (63) 

From  here  we  see  that  if  a=p  or 

9  =  1 — L*  0.292893 186  (64) 

then  the  method  is  second-order  accurate.  If,  however,  neither  of  these  conditions  hold 
then  the  method  is  only  first-order  accurate.  Thus,  a  reasonable  tradeoff  between  stability 
accuracy  seems  to  be  the  choice  given  by  Eq.  (64).  Also,  it  is  convenient  to  pick  a  and  p 
so  that  the  same  matrix  is  obtained  for  each  integration  step  in  (58).  This  means  that 
<x£=p (f-2fl)  from  which  we  may  determine  (requiring  a>p) 

a=__  and  0  =  — .  (65) 

For  the  value  of  0  given  by  Eq.  (64)  the  values  of  a  and  p  are  then 

a=  2 -V2  *0.5857864380  and  /?  =  V2  - 1 »  0.4142135620.  (66) 

This  gives  p/a=1/V2»0.707107. 

3.4  Numerical  Integration  of  the  Navier-Stokes  System 

The  application  of  the  9-scheme  for  time  integrating  the  discrete  Navier-Stokes  system  is 
outlined  here.  For  a  more  complete  description  see  reference  [13].  The  discrete-time 
values  for  the  dependent  variables  are,  as  before,  denoted  with  a  subscript.  For  example, 
the  x-directed  flow  is  given  by 


u„°=u°(nAr).  (67) 

The  other  dependent  variables  are  time  indexed  similarly.  We  assume  that  the  values  of 
the  dependent  variables  are  given  for  n=0,  that  is,  we  have  a  divergence-free  flow  from 
which  to  start.  Again  introducing  the  numerical  parameters  9e(0,1/2)  and  ct,pe(0,1)  with 
a+p=1 ,  we  must  solve  the  following  three  systems: 
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[4]U"v;— +  avtgtf] u°*+e  -UtflpUe  -[E^K+e  = 
u  At 

f„l+e  +  w*(u«+<»v*+e)  -M^]u®  -([C$]u„)u„  -([C|t]v„)u„ 

[4-]V"nf—  +  qvt^]v^ -[Djj^Le  -C^K+e  = 

C7  Af 

Co  +  Wr(uB+0,  vB+0)  -  M^]v®+0  - ([C|]u„)vn  - ([C^]v„)v„ 
[^]Tu»+e+[4]Tvn°+0=O  (68) 

[ Ey  ]  un+0  =  cn+e 

[^]Tv“+e  =  dn+e 

U«+e  =  Uinf 

v®+e  =  vinf 

for  u^e"0,  v^e°,  pm9,  ln+o,  and  iw  next  solve 

[ ^'j ]  "(1-26)  Ar+e  +  ^V[^<y]Un+l-9  +  ([^-^n+l-e^n+l-e  +  ([Qi]Vn+l-0)U;i+l-0  = 

fn+i-0  +  W  (u„+[_g,vn+|_^)  —  av[5(y]un+g  —  [Djj  ]pn+g  —  [•£</• ]^*n+0 

[Ay]  V(i':'2g)VA7'  +  ^v[^.y3v”+i-8  +([C^]u„+1_0)v„+1.0  +([C|t]vn+1.,)v„+I_e  =  (69) 

f*+i-<9  +  WK(u®+i_0,  v®+1_0)  —  av[^]v„+0  —[A/  ]P>j+0 

B 

Urt+l-0  —  ^inf 
vn+!-0  “*  vinf 

for  u^-e0  and  v^14°;  finally  solve 

-p/]pi+,  -[^]Xn+1  = 

f„+l  +  W^(u^+1,VB+1)  —  ^v[5^]u®+1_e  —  ([C,j/t]un+l-e)u/i+l-fl  —  ([Q/*]vn+l-e)Un+!-9 

[4]V°+'~^;+1-+av[%]v^1  -[^ri*i+i  -w«+>  = 

fn+i  +  Wr(uB+1, vB+i)  -  ^v[5^]v®+|_e  -  ([Cyi]un+i-0 )vn+]_0  —  ([Cy*]vn+i_e)vn+i-e 
[Dy  ]Tu®+i  +  [Dj]  ]Tv®+1  =  0  (70) 

[^•]T«®+i=c„+I 

[£*]TvL  =dn+1 

^n+l  — 

^n+l  —  ^inf 
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for  u^i°,  v^°,  pm1,  A.n+1,  and  (w  For  the  particular  choice  of  integration  parameters 
given  by  Eqs.  (65)  and  (66)  (i.e.,  9=1  -1/V2,  a=2-V2,  and  p=V2-1),  the  e-scheme  for  the 
Navier-Stokes  system  seems  unconditionally  stable  [14]. 


3.4.1  Solution  of  the  Stokes  Problem 

For  each  time  step,  two  Stokes-type  systems  and  one  advection-type  problem  must  be 
solved.  The  Stokes  systems  include  the  incompressibility  condition  which,  in  the  discrete 
system,  is  an  algebraic  constraint.  Here  we  develop  an  iterative  method  for  solving  the 
Stokes’  systems;  it  is  based  on  a  one-shot,  conjugate-gradient  algorithm  with 
preconditioning.  The  method  is  driven  by  the  Lagrange  multipliers  p„,  X„,  and  [15].  The 
solution  of  the  advection  subproblem  is  outlined  in  the  following  section.  We  begin  by 
describing  the  direct  solution  to  the  Stokes  problem  and  its  relationship  to  some  simple 
iterative  schemes. 

The  discrete  Stokes  system  may  be  written  in  the  form 

*  ’  5[AiJ]u  +  v[BiJ]u  =  fx+[Di?]t>  +  [Eijfr 

8[AV ]v  +  v[Bv ]v  =  fy  +  [Dy  ]p  +  [Ev ]p 
[Dy  ]Tu  +  [£),J ]Tv  =  0  (71) 

[^.]Tu  =  c 
[£*]Tv  =  d 


where  8  is  the  reciprocal  of  the  time  step  and  where  we  have  generalized  f  and  fy  to 
include  the  entire  right-hand  sides  of  Eqs.  (68)  and  (70).  This  is  the  form  on  which  we 
concentrate. 

In  Phase  I  we  used  direct  methods  to  solve  the  Stokes  problem.  Since  we  were 
considering  small  systems  for  the  purpose  of  development,  this  presented  little  problem. 
However,  when  considering  real-world  situations  the  system  size  will  typically  make  direct 
solution  methods  so  computationally  intense  as  to  place  a  bottleneck  on  real-time 
computation.  For  large  systems  it  is  more  practical  to  use  iterative  solution  techniques. 
Such  techniques  are  usually  based  on  the  principle  of  minimizing  a  functional  of  the 
dependent  variables.  The  process  is  designed  so  that  the  functional’s  minimum 
corresponds  to  the  solution  of  the  original  system.  Therefore,  in  this  subsection  we 
present  the  direct  solution,  a  general  iterative  scheme,  and  a  preconditioned  conjugate 
gradient  scheme  developed  by  Glowinski,  et  al  [14].  With  preconditioning,  he  reports  an 
effective  decoupling  between  system  size  and  number  of  iterations  until  convergence. 

3.4.1 .1  The  Direct  Solution 

The  direct  solution  to  (71)  may  be  found  by  forming  an  augmented  linear  system  with 
dependent  variable  (u,v,p,X.,n)T.  First,  to  simplify  matters  we  define  the  matrix  m 
according  to 


[Tg]mS[Ag]  +  V[Bg].  (72) 
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Note  that,  in  general,  [7"g  is  invertible.  Now  we  collect  all  the  equations  in  (50)  into  a  single 
matrix-vector  equation. 


'  [Ty] 

0 

~[Dy] 

-[Ey] 

0  ' 

fU] 

(f'] 

0 
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~[Dy] 
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~[Ey] 
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~lDy  ? 

-{Djjf 

0 

0 

0 
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0 

-[Ey? 
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X 

-c 

0 

V 

-[Ey? 

0 

0 

0  j 

Obviously,  the  direction  solution  of  the  above  would  not  be  the  most  efficient  way  to  solve 
the  problem  on  a  computer.  The  matrix  in  (73)  is  very  sparse;  indeed,  the  submatrices 
themselves  are  sparse.  Therefore,  when  solving  the  system  directly  it  is  wise  to  exploit 
the  structure  of  the  above  system.  However,  from  immediate  inspection  of  Eq.  (73)  we  do 
get  some  clues  as  to  the  behavior  of  our  system.  For  one,  it  is  symmetric,  or  self-adjoint. 
Therefore,  the  eigenvalues  will  be  real  and  distinct,  assuming  the  system  is  of  full  rank. 
We  knovy  from  the  governing  equations  that  the  pressure  p  is  only  defined  to  within  an 
arbitrary  constant,  so  it  might  be  the  case  that  the  above  is  not  of  full  rank.  However,  we 
show  below  that  the  matrices  [D*],  [D?],  and  [Eg  are  of  full  rank  and  therefore  (73)  has  a 
unique  solution. 

Since  the  matrix  [7"g  is  invertible,  we  can  solve  u  and  v  in  terms  of  the  remaining 
dependent  variables  using  the  first  two  rows  in  (73).  The  result  is 

u  =  [TyTlfx  HTyT'W^  +  iTyT^X, 
v  =  [^r'fy  +[Ty]-l[D?}p  +  [Tyr'[Eij]ix.  (,H) 

By  substituting  these  values  into  the  forth  and  fifth  row  equations  we  get  the  equations  for 
the  control  constraints 


[EgflTyT1  f*  HE9f[T9r\Dfa  +  [Ev]T[T„r'[E„y.  =  c, 

[Ey  ]r[TyT'fy +[Ey]T[Tv  ]"*  [D*  ]p  +  [Ey  ?  [Ty  ]-'  [Eg  ]p  =  d . 

The  matrix  [Eg  is  of  full  rank  since  for  every  /'we  can  find  a /such  that  [HiFflOO  (Le.,  the 
boundary  rc  intersects  the  domain  Q  somewhere).  Therefore  the  matrix  [Eg‘[Eg  is 
square,  of  full  rank,  and,  thus,  invertible.  Likewise  since  [Tg'1  is  of  full  rank,  the  matrix 
[EgT[Tg'1[Eg  is  also  invertible.  Using  a  similar  arguments  we  can  make  the  same 
conclusions  about  the  matrices  [D*[[T$\D*\  and  [D??[T^[D?].  This  motivates  the 
convenient  definitions 


[Zy}  =  [Ey}T[Ty]-'[Eg], 

[&•*]  =  [Dy  ]T [TyYX[E>y  ]. 

[A  Yy]m[Dry\T[TyY\Dry] 

[A^.]S[Af]  +  [Ay]. 

Using  the  fact  that  all  the  above  matrices  are  invertible,  we  may  begin  solving  for  the 
dependent  variables  in  (73).  From  Eq.  (75)  we  may  solve  for  X  and  p  in  terms  of  p 
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n = -[%  r  [Ey  fiTf  r'lDji  ]p + [%  ]_1  d  -  [s^,.  ]_1  ]rt^  rl  f y . 


Finally,  the  pressure  p  is  found  from  the  third  row  of  Eq.  (73)  and  the  above  relations.  We 
have 


P  =  ]r[7// r* [Eylfa- +  [Dy  ]r[7),- ]-1  [Etj ]n  +  [Dy  ]r[7)J- ]“'  f *  +  [D-  ]T [Ty ]“' f y ).  (78) 

Equations  (74),  (77),  and  (78)  may  be  back  substituted  to  find  independent  expressions  for 
all  the  dependent  variables. 

3.4. 1 .2  Overview  of  Iterative  Algorithms  for  Linear  Systems 


Typically,  iterative  schemes  for  solving  linear  systems  can  be  formulated  in  terms  of  the 
minimization  of  a  functional.  For  example,  consider  the  simple  matrix-vector  system 

[4ylx  =  y  (79) 


Where  x  is  the  unknown  vector  and  y  is  given.  We  shall  assume  that  [4j  is  symmetric 
and  positive  definite.  The  solution  is,  of  course,  x=[A;j-'y.  However,  rather  than 
computing  this  directly  we  choose  to  minimize  the  functional 


/(x)  =  ^xr[4]x-xry  +  c,  (80) 

where  c  is  some  constant  vector.  Notice  that  since  [Ail  is  positive  definite  the  solution  to 
the  minimization  problem  is  x=[Aj]~1y,  which  is  found  by  direct  differentiation.  The  idea 
here  is  to  iteratively  update  the  current  solution  vector  xm  according  to  some  (intelligent) 
process  guided  by  the  minimization  of  f{%).  More  specifically,  we  update  xm  according  to 
the  formula 


xm+ 1  xm  (81) 

where  dm  is  the  current  search  direction  and  a^>0  is  the  search  length  in  the  direction  dm. 
The  art  of  the  algorithm  is  in  the  choice  of  the  search  direction  dm.  The  search  direction 
may  be  found  by  direct  scalar  minimization  of  Eq.  (80).  Specifically,  define  the  scalar 
function  cp(am)  according  to 


<i>(am)  =  f(Xm  -“A) 

=  —n2iiT 


=  ja2mdTm[AiJ]dm -amdrm[AiJ}xm+amdlmy  +  ^Tm[Aij]xm  -x‘my  +  c 


(82) 


where  we  have  used  the  fact  that  [AJ  is  symmetric  (as  in  the  fluids  case).  By 
differentiating  <p(am)  with  respect  to  am  and  setting  the  result  equal  to  zero,  we  find  the 
minimizing  value  of  amto  be 


dm([A]x"~y) 

<[Ay]dm 


(83) 
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From  the  above  we  see  that  the  value  of  am  becomes  zero  once  the  xm=x  is  reached. 
This  has  an  interesting  interpretation  in  the  inner  product  space  generated  by  [Ail-  Define 
the  inner  product  (•,•)*  and  induced  norm  I  •  L  by 


(*>y)A  =xr[4/]y, 

and 


kL  =(x,x)A 


1/2  _ 


(84) 


Since  we  have  assumed  that  [A,j]  is  symmetric  and  positive  definite  we  know  that  this  does 
indeed  form  a  valid  inner  product.  Now  consider  the  unit  vector  in  the  direction  of  dm  given 
by 


d„  = 


11/2  • 

I  [  Ay  ]dm 


(85) 


-  » 


In  this  context,  a^dm  may  be  expressed 

amdm  —  dm  —  y) dm  —  (dm  ,xm  —  [Ay]  v)  d„ 


(86) 


so  that 


*m+ 1 


~xm  _(dm  ,Xm  -[Ay]  'y)dm  » 

=  xffl  +  Proj^  ([Ay  r1  y  -  Xm  )dm  . 


(87) 


Thus,  we  see  that  at  each  iteration  we  are  adding  in  the  offset  of  xm  (from  the  solution 
value  x=[Af ’y)  along  the  search  direction  dm.  Consider  if  the  search  directions  {dm}  could 
be  chosen  so  that  they  were  orthogonal  in  the  inner  product  <•,•)*  We  say  that  the  {dm}  are 
[/^-orthogonal  or  “[^-conjugate”.  Then  {d^"  would  form  a  basis  of  the  inner  product 
space  and  the  algorithm  would  be  guaranteed  to  converge  in  N  steps,  where  N=size(xm). 
This  tenet  is  the  foundation  of  the  conjugate-gradient  algorithms. 

3.4.1 .3  A  Functional  Equation  for  the  Stokes  System 

Now  consider  the  discretized  Navier-Stokes  system.  Notice  that  if  we  know  (p,X.,p.)  then 
we  can  find  (u,v)  by  direct  substitution  from  Eqs.  (53).  This  suggests  an  iterative  scheme 
for  solving  the  above  system.  We  can  create  a  functional,  which  when  minimized,  yields 
the  solution  to  the  linear  system  of  Eq.  (73).  This  functional,  f,  will  be  a  function  of  the 
variables  (pA,n)  only.  To  exemplify  this  idea  we  simplify  the  situation  by  removing  the 
solid  object  in  the  fluid  domain  so  that  A,=p=0.  In  this  case  we  find  from  Eq.  (78)  that 

p  =  -([Af]+[A $\Df]T[T„Tl1x  HD^fiTyT1  fy).  (88) 

However,  we  wish  to  find  this  solution  iteratively.  So,  in  the  spirit  of  Eq.  (80),  we  form  the 
functional 
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[Tij]  0 
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T 

ff1 

/(u,v,p)  =  - 

V 

0  [Tjj\  -[D*  ] 

;:=:VV 
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fv 

z 

IP, 

-[Aff  -P|]r  o  , 

3 

t°J 

=  iur[^]u  + 1  vr[7:.]v  - ur[Z)/]p  -  v^Djjp  - urf x  -  vrf  y  . 

We  assume  that  the  augmented  matrix  is  positive  definite  so  that  a  minimization  of  /(u,v,p) 
has  a  unique  solution.  (This  condition  relies  upon  the  relative  magnitudes  of  8  and  v.  As 
long  as  8  is  sufficiently  large,  indicating  a  small  time  step,  then  this  will  be  the  case.)  It  was 
mentioned  previously  that  we  could  substitute  out  u  and  v  using  Eqs.  (74)  making  f  solely 
a  function  of  p,  however,  to  avoid  the  algebraic  clutter  we  shall  refrain  from  doing  so  at 
present.  Just  remember  that  u  and  v  are  implicit  functions  of  p.  Consider  the  following 
iterative  updating  scheme  for  pressure  pm 

Pm+1  =Pm-«mdm  (90) 

-  f 

where,  as  before,  dm  is  the  current  search  direction  and  am  is  the  current  search  length. 
From  Eqs.  (73)  (with  X,=p=0)  we  find 

um+l  =[^r1fX+[^r'P/]pm+.  =lTij]-'fx+[TyT\D,f]pm-am[TyTl[D1f}dm> 

Vm+1  =[^r‘fy  HTyT'lD?] pm+,  =[Tvr'fx  +[Tur'[D?]pm -aJ.TgY^lD^. 

Equations  (90)  and  (91)  constitute  the  update  scheme  for  an  iterative  algorithm  once  we 
have  decided  upon  the  search  directions  dm.  The  search  length  am  can  be  found,  as 
before,  by  direction  minimization  of  Eq.  (89)  with  respect  to  am.  To  do  so  note  that  Pm<-i, 
Umn,  and  v,„m  are  all  functions  of  am  through  Eqs.  (90)  and  (91).  Their  derivatives  with 
respect  to  am  are  given  by 


<?Pm+l 


=  -d„ 


dam 


=  ~[Ty  ]-'[£)/  ]dm, 
=  -[Ty]~'[D*]dm. 


(92) 


Defining  the  scalar  function  cp(am)  as 

^(®m)  sy^(Pm+l(®'m)’**m+l(®m)’'rm+l(®m))  (9^) 

we  take  its  derivative  and  set  it  equal  to  zero 

d<p(am)  1  ^(pnn-l’unn.|’viin-l)  ^m+l  ^  ^(Pm+l  ’  uwtl » vn»-l  ))  d  vm+l  (94) 

do.m  &Pm+l  dtlm+i  d&m  ^vm+l  ^am 


The  above  equation  is  solved  for  am  with  the  result 
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(95) 


am  =' 


^[Ag3pw+dl([D^]r[rffr1f1+[gj3r[rffrlfy) 

d„[A y]dm 


Note  that  this  formula  is  not  as  convenient  as  the  one  generated  for  the  simple  system 
[AJx=y.  To  use  the  above  equation  we  must  invert  and  store  the  matrix  [Tjj  and  form  the 
matrix  [Aj,  a  costly  procedure  considering  the  potential  size  of  the  system.  However,  the 
procedure  would  have  to  be  done  only  once  for  each  value  of  the  parameter  8. 
Depending  upon  circumstance  it  might  be  faster  simply  to  solve  for  the  pressure  directly 
using  Eq.  (88). 

A  simple  interpretation  results  if  the  system  is  such  that  the  drive  vectors  f  and  fy  are 
representable  from  a  single  vector  g  in  the  pressure  domain  according  to 


-  * 


fx  =[Z>/]g, 
fy=[^]g- 


(96) 


Then  the  update  formula  for  the  pressure  p^i  simplifies  to 


Pm+I  Pm  (dm>Pm  "*"S)Adm  (97) 


where  the  notation  is  analogous  to  that  described  in  (84).  In  this  case  the  solution  for 
pressure  is  p=-g  which  comes  directly  from  Eq.  (78)  and  is  implied  by  Eq.  (97). 

3.4.1 .4  A  Simple  Iterative  Algorithm  for  the  Stokes  System 

Here  we  present  a  iterative  algorithm  which  will  solve  the  Navier-Stokes  system  for  the 
case  when  X=p=0.  This  is  a  simple  "ad  hoc"  method  where  the  results  of  the  previous 
subsection  are  essentially  ignored.  That  is,  the  choice  of  am  is  not  guided  by  any  direct 
test  of  the  dependent  variables.  Instead,  we  are  guided  to  this  algorithm  from  the  physics 
of  the  fluid-flow  problem.  This  presentation  shows  that  it  is  possible  to  implement  very 
simple  algorithms  based  on  physical  intuition  (that  are  guaranteed  to  converge),  however, 
the  rate  of  convergence  may  not  be  practical. 

This  method  chooses  the  search  directions  for  pm  as  the  divergence  of  the  velocity  field. 
That  is,  we  choose 


Am=[DhT  u, 


+  [DhTy 


(98) 


which  is  the  discrete  analogue  of  d=V  U.  The  values  of  jpm  um  and  vm  are  then  updated 
according  to 


Pm+l=Pm-adm  (") 


and 


um+,  =[T9rlF+[T„Yl[D?] pm+i, 

vm+,  =[^r1fy+[7),'r1[Aj]Pm+i. 
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which  are  the  same  as  Eqs.  (90)  and  (91)  except  for  the  search  length  parameter  a™. 
Here  all  the  search  lengths  am  are  chosen  to  be  the  same  number,  a,  which  lies  in  the 
interval  (0,2v/A/).  With  these  conditions  the  algorithm  is  guaranteed  to  converge  [1 6].  To 
see  that  the  convergent  solution  is  indeed  the  correct  solution  note  that  the  algorithm  is 
essentially  a  fixed-point  algorithm  which  converges  when  the  (discrete)  divergence  of  the 
flow  is  zero.  Thus,  at  that  time  the  search  direction  will  be  zero.  The  value  of  the  fixed 
point  (p,u,v)  is  then  found  from  Eqs.  (98)  to  (1 00) 


p  =  p-ad  =  p-  a([Z),f  ]ru  +  [Z>J]v) 

=  P  -a([Z>/  f[TyT'fx  +[Aj][^r'fy  +[A*]p) 


(101) 


whose  solution  is 

P=-[^r,([^f[2;yrlfx  (102) 

From  Eq!  (88)  this  is  seen  to  be  correct. 

The  above  algorithm  will  be  quite  slow  for  flow  at  large  Reynolds  number  (i.e.,  a  small 
value  of  viscosity  v).  To  further  aggravate  this  situation,  the  time  step  At  must  be  small  to 
follow  the  fast  dynamics  of  these  flows.  Thus,  to  improve  the  speed  of  convergence  we 
use  a  preconditioned  conjugate  gradient  version  of  the  above  algorithm. 

3.4.1 .5  A  Preconditioned  Conjugate-Gradient  Algorithm 

Here  we  outline  a  conjugate-gradient  algorithm  with  preconditioning  which  has  been 
presented  by  Glowinski  et.  al.  [14],  The  method  is  initialized  by  picking  (arbitrarily)  vectors 
p0,  X0,  and  po.  Then  we  solve  the  linear  system 


8[Ay]  u0  +  v[Bv]  u0  =  f*  +[£>/ ]p0  +[£*  ]*-o 
<5[4-]v0  +  v[5,y]v0  =  fy  +  [£>J  ]p0  +  [£(>]Po 


for  Uo  and  Vo.  Now  compute  the  initial  "remainder"  vectors 

s0  =[£>i?  ]Tu0+[Ay]Tvo> 

rg  =[£,y]Tu-c,  (104) 

ry0  =[^]Tv-d. 


To  continue  we  must  solve  two  auxiliary  systems  which  may  be  identified  as  the  ''costate" 
equations  for  pressure  and  for  the  fictitious  domain  Lagrange  multipliers.  Beginning  with 
the  definitions 

♦  <t>'  =(<t>i\ieJx  ),  (105) 
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the  discrete  system 


[By  ]<|>o  ~[Ny  ]<t>o  =  s  -  [By 
<t>oS  =° 


(106) 


where 


8Gt  dGj  8Gj  dGf 
dx  8x  +  dy  dy 

dQ. 

for  i,j  eJl 

b‘‘L 

~  8  Gt  dGj  8Gt  dGj~ 
dx  dx  +  dy  dy 

dQ 

for  i  eJuj  eJB 

(107) 

£ 

II 

J“1  —  * 

‘ 8Gj  8Gj 

17 n* +  dy  "'J 

GidT 

for  i,j  eJ\ 

constitutes  the  pressure  costate  and  must  be  solved  for  <po.  Note  that  cpo  is  solved  in  the 
pressure  space  P"(! Q ).  The  above  system  is  the  discretization  of  the  weak  form  of  the 
Poisson  equation 

-V20  =  V-U  inQ 

^  =  0  onT0  (108) 

8n 

0  =  0  on  r, 

which  typically  occurs  somewhere  in  the  formulation  of  incompressible  fluid  simulations. 
The  second  auxiliary  system  to  be  solved  is  the  costate  system  for  the  Lagrange 
multipliers  X  and  p.  In  this  system  we  solve 

[Mfj}g 5-r0x  and  [Mf] |g$ (109) 

for  the  vectors  go"'  and  g0y.  System  (1 09)  is  the  discrete  version  of  the  equation 

a(g,b)  =  j*  r-bdT  Vbe5(rc)  (110) 

where  a(v)  is  some  bilinear  form  on  B( rc).  We  shall  choose  a(v)  to  be  the  standard  inner 
product  so  that  the  matrices  [M*\  and  [MilY]  are  simply  the  identity. 

From  the  solutions  to  the  two  auxiliary  systems  we  form  the  initial  "gradient"  vector 
go=(gop,gox,goy)  (corresponding  to  p0,  Xq,  and  po)  according  to 


'to' 

r<5  0+^0 1 

OTQ 

o 

III 

gSo 

w 

■ 

✓ - 

*i  jrt 

°Vc  °  * 

___ 

Finally,  we  set  our  initial  (negative)  search  directions  w0=(wop,w0x,Woy)  equal  to  the 
gradient  vectors  (gop, gox-9oy)  or 
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(  n\ 

"go  ^ 

W-J, 

ii 

g» 

w 

w 

(112) 


This  completes  the  initialization  of  the  algorithm. 

For  the  main  loop  of  the  algorithm  we  assume  that  the  iteration  step  m> 0  and  that  pm,  K,, 
Mm.  um,  vm,  gm,  and  wm  are  all  known.  To  find  the  next  solution  iterates  pm+i,  Vn,  Mm+i, 
Um+ii  vm+1,  gm+i,  and  wm+i  we  proceed  as  follows: 

Solve  the  linear  system 

S[Ay  ]u  +  V[By  ]u  =  [D*  ]w£  +  [Ey  ]< 

5[4.]v  +  v[^]v  =  [DjK  +[Ey]<  {]lh 
for  theinferim  velocities  u  and  v.  Now  form  the  interim  remainder  vectors  s,  r*,  and 

s  =  [D/]tu  +  [D>:]tv, 

rx  =[£■,■,■  ]Tu,  (114) 

ry=[Ey]Tv. 


and  set  the  gradients 

g£=<5*  +  v«,  g*  =  rx,  g:=r\ 
where  cp  solves  the  system 

<j>s  =  0 

(the  definitions  for  cp1  and  cp8  are  as  before).  Next  compute  the  search  length  parameter 
pm  according  to 

(Or[4f]g£  +(r*)rU/]g£  +«)r[4y 

Q  —  - — — - - - -  (li/J 

where  the  matrix  entries  are  given  by 

Ay  =  j^GjGj  dCl  for  i,j  eJ  and  Ay=^HiHJdT  for  i,j  eK.  (118) 

Update  the  vectors  using  the  search  directions,  previous  iterates,  and  search  length 
parameter. 


(115) 


(116) 
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P  m+1  P  m  Pm  ^  m 


«+l  =  ^~PmWm 
M’m+l  “  Pm  “  PmWm 

(119) 

Um+1  Um  PwU/n 

^m+l 

Sm+!  —  Sm~"Pm8m 


Finally,  if 

(rm+l)r[-'4#  ]§m+l  +(rm+l)  [ -^y  ]8m+l  +  (rm+l)  My  ]§m+l  (120) 

(roP ) T[Ay  ]§o  +  (ro  )TIAj  ]§o  +  (ro  ) T[Aj  l6o 

where  s  is  some  solution  tolerance  parameter,  then  stop  the  algorithm  and  take  the  last 
iterate^  the  solution  to  the  Stokes  problem.  Otherwise,  compute  the  search  direction 
parameter 


7m 


)  [Ay  ]gm+\  ■*"(rm+[)  [Ay  ]§m+l  +  (rm+l  )  My  ]§m+I 
{<)T[Ap]gpm  +K)rlAif]gxm  +  (r>)T[Aif]gym 


(121) 


then  update  the  search  direction  wmto 

wm+1  =gm+ymwm.  (122) 

Set  the  iteration  counter  to  m=rm  1  then  restart  the  algorithm  from  Eq  (113). 

3.4.2  Solution  of  the  Advection  Problem 

The  discrete  advection  problem  may  be  written 

8[AV  ]u  +  v[B0-  ]u  +  ([C&  ]u)u  +  ([Cy  ]v)u  =  fx  ^ 

8[Ay  ]v  +  v[B;j  ]v  +  ([<$  ]u)v  +  ([C£  ]v)v  =  f  y 

where,  again,  5  is  the  reciprocal  of  the  time  step  and  f*  and  fy  have  been  generalized  to 
include  the  entire  right-hand  side  of  Eq.  (69).  This  is  a  nonlinear  system  due  to  the  terms 
involving  the  tensors  [C,/]  and  [C/],  however,  there  are  no  constraints  and  the  pressure 
p  and  Lagrange  multipliers  X  and  p  do  not  appear.  This  is  set  of  second-degree  algebraic 
equations  in  u  and  v  where  there  are  as  many  equations  as  there  are  unknowns. 
Therefore,  we  can  use  standard  techniques  for  solving  this  nonlinear  set  of  equations.  In 
particular  Newton’s  method  is  readily  applicable.  Yet,  it  is  also  possible  to  use  other 
iterative  techniques  such  as  a  nonlinear  conjugate  gradient  method.  To  do  so,  however,  a 
least-squares  formulation  of  Eqs.  (123)  must  be  provided. 

3.4.2.1  A  Functional  Equation  for  the  Advection  System 

A  functional  equation  for  the  advection  problem  may  be  written  just  as  was  done  for  the 
Stokes  problem.  Presumably  then,  minimizing  the  functional  with  respect  to  u  and  v  will 
result  in  the  solution  to  the  advection  problem.  However,  in  this  case  the  functional  is  of 
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third  order  and  not  necessarily  convex.  Thus,  the  minimization  will  require  more 
sophisticated  techniques.  Proceeding  analogously  as  in  Section  3.4.1 .3,  the  resulting 
functional  f  is 


/(U,  V)  =  u  t[T0  ]u  +  v  r[Tij  ]v  +  u 


T  ([C^.]u)u+([C|t]v)]u  +  v7'  ([C|]u)u+([C^.]v)]v-urfx  -vrfy 


(124) 


where  the  matrix  [TJ  is  as  before  (defined  by  Eq.  (72)).  As  of  yet  we  have  not  pursued  this 
solution  strategy.  Instead,  we  have  relied  upon  a  linearized  approximation  of  the 
advection  system.  It  may  be  necessary  to  solve  the  full  nonlinear  system  if  we  encounter 
any  accuracy  problems  in  the  future. 

3.4.2.2  A  Linearized  Approximation 

Currently,  we  solve  this  problem  using  a  linearized  approximation  to  the  second-degree 
terms  in  the  advection  equations.  We  do  this  by  using  the  previous  iterates  of  u  and  v  for 
the  initiahmultiplications  of  the  advection  tensors,  that  is 

([Q*  lUn+i-e^n+i-e  +  ([C|t ]v«+i-e)u®+i-0  «([C,^]u®+e)u„+i_e  +([C^ ]v„+0)un+1_9  , 

and  (125) 

]u/i+l-0  )V«+l-0  +  ([  cfjk  ]Vn+l-tf)V;i*l-0  ~  ]U«+0  )vn+l-0  +  ^ijk  ]vn+0  )v/i+l-0  • 


It  has  been  noted  that  there  is  practically  no  loss  in  accuracy  and  stability  when  using  this 
approximation  [13].  The  resulting  linear  system  is  substantially  smaller  than  the  Stokes 
system.  Consequently,  using  direct  techniques  to  solve  this  system  has  not  produced  any 
significant  bottleneck  concerning  real-time  computation. 


4  Control  Problem  Structure  for  Fluid  Flow 


We  have  been  investigating  the  control  problem  structure  for  fluid  flow  in  both  the  discrete 
(FEM/FDM)  and  the  continuous  case.  In  the  discrete  case  our  system  is  governed  by  a 
set  of  nonlinear,  first-order,  ordinary  differential  equations  with  constraints.  We  consider 
aspects  and  techniques  of  controlling  these  equations  directly.  However,  this  approach 
does  not  necessarily  provide  us  with  a  clear  technique  for  controlling  the  actual  physical 
system.  Thus,  we  are  also  investigating  the  boundary  layer  control  of  the  continuous, 
incompressible  Navier-Stokes  equations. 

4.1  The  Discrete  Case 

The  control  theory  for  the  discrete  Navier-Stokes  equations  is  in  the  very  earliest  stages  of 
development.  Our  approach  has  been  influenced  by  the  important  work  of  Glowinski, 
Temam,  Lions,  Banks,  and  others.  However,  for  the  most  part  this  work  has  not  led  to 
“practical”  feedback  laws  that  can  be  tested  in  simulations  or  in  the  laboratory.  Since  this 
is  our  primary  goal,  our  approach  has  accordingly  been  somewhat  pragmatic -find  control 
laws  based  on  measurable  quantities  that  can  be  implemented  with  reasonable  actuation 
schemes  and  evaluated  in  simulation.  We  have  limited  ourselves  to  finite  dimensional 
measurements,  and  to  actuators  with  finite  dimensional  influence  functions.  Thus,  we 
assume  that  the  flow  state  can  be  measured  only  at  a  finite  number  of  points  (near  or  on 
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the  boundary  of  the  object  in  the  flow).  And  we  assume  that  the  actuators  act  at  a  finite 
number  of  points. 

To  simplify  the  presentation,  consider  the  numerical  representation  of  the  system 
dynamics  in  the  following  form: 

[AJJ  +  (t[5]]C/)C/  +  [C]U  + 1 D]P  =  F  +  [Q]  A 
[DT]U  =  0 
[Qt]U  =  X 

where  [A],  [C\,  [D\  and  [Qj  are  matrices  and  [[fij]  is  a  rank-three  tensor.  The  quantity  U  is 
the  discrete  vector  representing  the  velocity  field,  the  vector  P  is  the  pressure,  the  vector  A 
is  the  fictitious  domain  Lagrange  multiplier  for  an  embedded  object,  the  vector  F  is  the 
applied  force,  and  vector  X  represent  the  boundary  control.  Since  the  velocity  and 
pressure  grids  may  be  different,  the  dimensions  of  U  and  P  are  not  necessarily  related. 
The  first  equation  is  the  discrete  approximation  of  the  momentum  equations.  The  last  two 
equations  represent  the  divergence  condition  and  the  boundary  control  condition  on  the 
object.  In  effect,  this  is  a  kind  of  differential-algebraic  system  -  a  control  system  with  an 
algebraic  constraint.  There  are  several  ways  to  handle  such  systems.  Our  approach  is 
based  on  singular  perturbation  methods  (suggested  by  the  numerical  methods).  We 
introduce  additional  state  variables  and  re-write  the  equations  in  the  form 

[A]U  +  ([[fijFX7  +  \!°V  +  \P\P  =  F+ [2]A 

^=at;+[DTJJ 

et=fc-[QTV  +  X 

Here  (s,n)  are  small  (positive)  parameters,  (a,p)  are  matrices  having  eigenvalues  with 
negative  real  parts,  and  (q,Q  are  state  variables  of  appropriate  dimensions.  The  state  £ 
represents  the  “dynamics”  of  the  boundary  control  system  (as  a  dynamical  compensator), 
and  the  state  £  is  an  artifice  introduced  to  approximate  the  divergence  condition.1  The 
parameters  (e,n)  are  “small”  indicating  our  intent  to  use  (singular)  perturbation  methods; 
but  also  indicating  the  physical  interpretation  that  the  dynamics  of  the  boundary  control 
system  are  “fast”  relative  to  the  (large  scale)  dynamics  of  the  fluid  flow.  In  effect,  the  flow 
conditions  near  the  boundary  change  slowly  (in  time)  relative  to  the  capability  of  the  control 
system  actuator  dynamics.  With  this  formulation,  several  options  are  available  for  control. 
Suppose  we  want  to  control  the  flow  to  stabilize  the  velocity  field  U  in  the  presence  of 
disturbances  F  acting  in  the  flow.  In  this  case,  the  control  objective  is  to  choose  the 
(boundary)  control  values  X{t)  to  “reject”  the  disturbances  and  stabilize  the  flow.  In  the 
absence  of  the  nonlinear  term  in  the  above,  this  is  an  elementary  problem  in  linear  control 
theory.  Since  disturbance  rejection  and  adaptive  stabilization  methods  have  been 
available  for  some  time  in  nonlinear  control  theory,  tools  are  available  in  principle  to 
resolve  the  stabilization  control  problem  for  the  Navier-Stokes  equation. 

4J2  The  Continuous  Case 

In  this  situation  we  consider  the  continuous  Navier-Stokes  equations  directly  in  order  to 
develop  some  physical  insight  into  the  controls  problem.  Perturbation  methods  are 


1  The  divergence  condition  state  is  not  necessary,  we  can  use  differential — algebraic  methods:  it  is  used  here  to 
avoid  the  complexities  of  this  condition.  In  control  theory  parlance,  this  condition  (without  the  state)  amounts  to 
sliding  mode  constraint. 
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employed  to  derive  approximate  differential  systems  that  describe  the  behavior  of  the 
boundary  layer  around  objects  in  the  fluid  domain.  In  the  spirit  of  Prandtl  we  proceed 
using  a  boundary  layer  expansion  of  the  Navier-Stokes  equations,  however,  we  wish  to 
retain  the  time  dependent  terms.  In  order  to  do  so  we  find  that  a  multiple  time  scale 
expansion  may  be  used  to  account  for  changes  in  the  normal  fluid  velocity. 

Using  the  principles  of 
differential  geometry,  we 
have  written  a  Mathematica 
program  to  symbolically 
compute  the 

incompressible  Navier- 
Stokes  equations  in 
arbitrary  coordinate 

systems  (see  Appendix  A). 

For  example,  consider  the 
airfoil  depicted  in  Figure  9 
embedded  in  the  2D  fluid 
domain.  Construct  a 
coordinate  system  around 
it  so  that  x  is  the  tangential 
distance  along  its  surface 
and  y  is  the  distance  in 
direction  normal  to  the  surface,  as  shown  in  the  figure.  (Since  we  are  going  to  be 
interested  only  in  the  boundary  layer,  y  is  defined  in  the  layer  even  if  the  airfoil  is  not  locally 
convex.)  The  surface  of  the  airfoil  is  described  by  the  curve  y=0. 

By  applying  an  asymptotic  boundary  layer  expansion  Q=(Re)y  in  the  normal  coordinate 
(i.e.,  a  stretching  in  the  y  direction)  then  taking  the  limit  as  the  Reynolds  number  Re  goes 
to  infinity  we  get  the  following  equations  for  the  dynamics  of  the  boundary  layer: 


Stretched  coordinate  £ 


Figure  9:  Airfoil  with  stretched  local  coordinates 


d  u  du  du  dp  dl  u 

- +  u - +  v - +  —-  = - - 

dt  d  x  <?£  dx  8L T2 


dp 


1 


du  +  ^ v  _  o 
d  x  d£ 


k{x)u2 


(126) 


where  u  is  the  velocity  in  the  x direction,  v  is  the  velocity  in  the  C  direction  (y  direction),  p  is 
the  pressure,  and  k(x)  is  the  curvature  of  the  surface  at  x  (in  the  2D  case  it  is  simply  a 
scalar,  in  3D  it  is  a  rank-two  tensor).  Notice  that  these  are  the  Prandtl  limit  equations  for 
the  flat  plate  with  the  addition  of  the  term  accounting  for  surface  curvature.  The  boundary 
conditions  are  given  by 

u(x,0;t)  =  0  u(x,C;t)—£ZT>ue(x,  0) 

v(jc,0;  t)-X{t)  v(s,C;r)  ^  >ve(x,0)  (127) 

p(x&t)  >pe(x,0) 


where  X(\)  is  the  boundary  control  (which  may  be  nonzero  only  on  portions  of  the  surface), 
and  Ue,  ve,  and  pe  are  the  inviscid  solutions  to  the  problem  (i.e.,  when  f?e-*»)  in  the  (x,y) 
coordinates,  the  so-called  outer  solution.  We  are  currently  investigating  these  equations 
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using  asymptotic  and  perturbation  analysis  along  with  numerical  methods  to  develop 
control  strategies.  One  point  that  is  immediately  obvious  from  the  second  equation  in 
(126)  is  that  pressure  gradient  is  highly  sensitive  to  transverse  velocity.  Thus,  as  we  know 
empirically,  at  points  of  flow  separation  there  are  very  large  pressure  differentials.  For 
example,  for  flow  over  corners  k(x)  becomes  very  large  and,  since  u  must  remain  nonzero 
by  Newton’s  first  law  (causing  separation),  the  pressure  gradient  is  also  quite  large.  Note 
that  Eqs.  (126)  involve  the  tangential  velocity  u,  the  normal  velocity  v,  and  the  curvature  k. 
We  have  been  attempting  to  use  the  normal  velocity  at  the  surface  as  a  control  parameter. 
Another  possibility  is  to  use  the  curvature  k  as  the  control,  as  in  compliant  membranes  or 
flap  actuators. 

One  final  point,  if  we  introduce  a  fast  time-scale  variable  x={Re)t,  we  can  recover  the  v 
time  dependence  in  the  second  boundary  layer  equation 


dv  dp  _  1 


k(x)u2 


(128) 


This  is' ah  important  equation  when  considering  control  strategies.  It  can  be  shown  that 
du/di=0  (to  first  order)  so  that  u  may  be  considered  constant  in  the  above  equation.  Thus, 
changes  in  normal  velocity  v  and  pressure  p  can  happen  on  very  fast  time  and  spatial 
scales.  This  equation  can  account  for  boundary  layer  "bursting"  seen  in  experiment. 


4.2.1  Boundary  Layer  Analysis  (Singular  Perturbation) 


Here  we  demonstrate  the  validity  of  Eqs.  (126)  and  (127).  Rather  than  developing  these 
relations  using  differential  geometry  as  in  the  Mathematica  code,  we  use  a  more  direct 
boundary  layer  approach  to  clarify  the  ideas.  The  results  of  the  differential  analysis  will 
simply  add  an  additional  term  accounting  for  surface  curvature.  We  start  with  the  two- 
dimensional,  incompressible  Navier-Stokes  equations  in  cartesian  coordinates. 


8u  8u  op 
u  +  u  —  +  v —  +  — =  v 
dx  dy  dx 

dv  dv  dp 
v  +  u  —  h-v —  +  —  =v 
dx  dy  dy 


f  d2u  ^  d2u 

r  d2v  d2v' 
dx2  dy1 


du  dv 
dx  dy 


=  0 


(129) 


where  the  over  dot  indicates  differentiation  with  respect  to  time,  u  and  v  are  the  fluid 
velocities  in  the  x  and  y  direction,  respectively,  p  is  the  pressure,  and  v  is  the  kinematic 
viscosity  (equal  to  1/Re).  Of  course  there  will  be  boundary  conditions  that  accompany  the 
above  system;  these  conditions  will  depend  upon  the  fluid  domain  and  the  particular 
material  composing  the  boundaries.  We  can  represent  these  boundary  conditions 
symbolically  according  to  the  following: 


B,[u(x,y;t),  v{x,y,t),p{x,y\t)}= for  i  =  1,2, 


Na 


(130) 


where  the  are  the  boundary  condition  operators  and  NB  is  the  number  of  boundary 
conditions.  Because  of  the  time  derivatives  we  also  have  initial  conditions  associated  with 
the  time-dependent  Navier-Stokes  equations.  That  is,  we  require  an  initial  state  in  time  for 
the  fluid  that  can  be  represented  as 
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(131) 


u(x,y;t)  U  =«i(x,y), 
v(.v,7;OUo=vi(-1c.>?)> 
p(x,yU)\t^=Pi(x,y), 

where  14  vu  and  p,are  the  initial  states  of  the  field  variables. 


We  are  interested  in  flows  with 
large  Reynolds  number  Re, 
thus,  we  proceed  by 
considering  the  limiting  case 
where  v->0.  For  the  moment, 
consider  the  problem  of  a  flat 
plate  immersed  in  a  uniform 
flow.  The  plate  is  infinitesimally 
thin,  aligned  in  the  plane  y=0, 
and  has  no-slip  boundary 
condition^.  This  situation  is 
depicted  graphically  in  Figure 
10.  When  v=0  we  have  the 
Euler  equations  whose 
solution,  in  this  case,  is  simple 
uniform  flow  which  violates  the 

no-slip  boundary  conditions.  Prandtl  observed  that  in  order  for  a  limiting  solution  to  make 
sense  there  must  be  a  “boundary  layer”  around  the  plate  where  the  flow  changes  rapidly 
from  zero  to  its  uniform  flow  velocity  [6].  Thus,  we  expect  a  rapid  transition  in  u  and  v 
around  the  positions  y=0.  This  idea  is  further  supported  by  the  fact  that  the  viscosity  v 
multiplies  the  highest  order  derivatives  in  the  above  equations,  suggesting  that  these 
derivatives  must  approach  infinity  in  the  limit  v-»0  in  order  to  satisfy  the  boundary 
conditions. 


Figure  10:  Uniform  flow  around  a  flat  plate 


In  what  follows  we  are  going  to  compute  an  asymptotic  expansion  of  the  Navier-Stokes 
equations  in  terms  of  the  viscosity  v.  The  solutions  to  the  expanded  system  will 
approximate  fluid  flow  for  small  viscosities.  Due  to  the  boundary  layer  around  the  plate  (or 
any  object  immersed  in  fluid  flow)  it  is  necessary  to  expand  the  equations  in  two  separate 
regions.  One  expansion  is  valid  in  the  boundary  layer  and  is  referred  to  as  the  inner 
expansion.  The  other  expansion  is  valid  in  the  region  outside  the  boundary  layer  and  is 
referred  to  as  the  outer  expansion.  Once  the  solutions  in  the  opposing  regions  are  found, 
they  must  be  matched  together  to  produce  the  full  composite  expansion.  This  matching 
process  can  be  tricky  when  consider  higher  order  terms. 


We  continue  to  focus  on  the  flat  plate  problem,  so  that  we  assume  the  existence  of  a 
boundary  layer  around  the  points  y=Q.  We  consider  the  case  of  wall  curvature  later.  The 
boundary  conditions  for  a  flat  plate  in  uniform  flow  are 


u(x,  0;  r)  =  0 
v(jc,0;  0  =  0 


for  x  >  0 


\im  u(x,y;t)  =  U^ 

y— >oo 

u(x,y,t)  =  U„  {Qrx<0  l]mv{xy.t)  =  0  (132) 

v(jc,y;O  =  0 

lim  p(x,y;t)  =  Px> 

y->oo 


where  Ux  and  Px  are  the  free  stream  flow  velocity  and  pressure,  respectively. 
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,4.2. 1 . 1  The  Outer  Expansion 


In  the  outer  expansion  equations  the  field  variables  are  denoted  with  capital  letters,  thus, 
we  have  U,  V,  and  P  as  the  outer  solutions  for  x-velocity,  /-velocity,  and  pressure, 
respectively.  To  begin  with  assume  that  the  field  variables  can  be  expanded  in  a  power 
series  in  v 


U(x,y:t)  =  U0(x,y;t)  +  vUl(x,y;t)  +  O(v2) 

V(x,y;t)  =  V0(x,y;t)  +  vV{(x,y,t)  +  0(v2)  (133) 

P(x,  y;t)  =  P0(x,  y;  t)  +  vPx  (x,  y,  t)  +  0(v2) 

Substituting  these  expansions  into  the  Navier-Stokes  equations  and  collecting  terms  with 
like  powers  of  v  yields  separate  systems  of  equations  and  boundary  conditions  for  each 
order  of  v.  The  zero-order  system,  0(1 ),  is  given  by 


u0+u0^  +  v0 

OX 


0(1) 


dx 


dU0 

A 

=  0 

dy 

dx 

■  dV0 

+  S n 

=  0 

0 

dy 

du0 

+  dVL 

=  0 

dx 

dy 

(134) 


0(1)  Bi[u0(x,y;t),V0(x,y,t),P0(x,y;t)]=l}l  for  condition  at  /->  oo  (135) 


U(l(x,y;t)\,^=ul(x,y), 

0(1)  V0(x,y;t)  Uo=v,.(x,/),  (136) 

P0(x,y,t)  Uo =Pi(x,y), 


while  the  first-order  system,  0(v),  is  given  by 


0(v) 


0,+u.&.+o,2!L+rM*r, 

dx  dx  dy  dy  dx  dx 


*.U.ZL  +  U&  +  V.2 yr,dr- 


dx 


dx 


dy 


d\ 


a2c 
dy2 
d\ 


dy  dy  dx  dy 


dUx  dV, 
dx  dy 


=  0 


(137) 


0(v) 

Bt [O, (x,  y;t),Vt (x, y,t),Px(x,y,t)]=  0  for  condition  at  /  oo 

(138) 

U\{x,y,t)  |,_o=0, 

0(v) 

vt(x,y,t)\l= 0=o» 

(139) 

^(x,/;r)  Uo=°> 
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The  zero-order  equations  are  the  Euler  equations  exactly.  Consequently,  it  is  common  to 
denote  the  field  variable  solutions  with  a  subscript  e.  That  is  our  variables  U0,  V0,  and  Po 
also  appear  as  14  14  and  p@,  respectively,  in  the  literature.  Note  that  once  the  zero-order 
(Euler)  system  is  solved  for  l/0,  V0,  and  P0)  the  first-order  system  is  linear  in  Ux,  Vh  and  Pu 

Since  the  expansion  of  (133)  must  hold  for  all  v,  the  initial  conditions  must  be  satisfied  by 
the  zero-order  system  while  the  first-order  system  will  have  zero  initial  conditions.  This 
fact  also  dictates  the  form  of  the  boundary  conditions  on  future  expanded  systems. 

4.2.1 .2  The  Inner  Expansion  (Boundary  Layer) 

We  shall  consistently  use  lower  case  letter  to  indicate  field  variables  in  the  inner 
expansion.  As  we  shall  see,  there  will  be  several  different  sets  of  variables  resulting  from 
the  application  of  multiple  transforms. 

In  boundary  layer  the  field  variables  experience  extreme  variation.  Thus,  we  expect  a 
rapid  variation  with  the  coordinate  y  normal  to  the  plate  around  the  points  y=0  (for  x>0). 
This  condition  suggest  the  stretching  transform  for  the  normal  coordinate  y 


dy~va  dQ 


(140) 


where  a  is  some  positive  number  to  be  determined.  By  converting  to  the  coordinate  C,  the 
boundary  layer  can  be  held  to  a  nonzero  thickness  as  the  viscosity  v  approaches  zero. 
Because  the  terms  of  the  Navier-Stokes  equations  that  contain  the  second-order 
derivative  o2/^2  have  a  factor  v,  the  proper  choice  for  a  is  known  to  be 

a  =—  (141) 

2 

If  it  is  taken  larger, we  ecounter  an  unbounded  limit,  if  taken  smaller  the  limit  equations  fails 
to  generate  any  new  information. 

The  divergence  condition  in  the  governing  equations  implies  the  existence  of  a  stream 
function  y /(x,y)  such  that 

v(,,j ,)— 042) 
dy  vl/2  dx  8x 

Since  we  expect  dy(x,y)/d& 0,  we  must  use  the  following  transformations  to  keep  u 
bounded  as  v-*0: 


¥'(x,^)s-^Ty/(x,vl/20, 

y 

u’(x,Q=u(x,v'nC), 
v'(*,0  =  ^77v(*,vl/2a 
p’(x,C)  =  p(x,vU20, 


(143) 
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where  the  last  three  transforms  follow  from  the  first.  Applying  the  stretching  transform  to 
the  Navier-Stokes  equations  along  with  the  above  transformations  for  the  field  variables 
yields  the  new  system 


du'  ,du'  dp'  d2u  d2u 
dx  dC,  dx  dx2  d£2 


.■  A 


,  dv'  ,  dv' 


ox 


dp'  2  av  av  ..... 

+  — =  1/ — r  +  v — T  (144) 

ac  a^2  ac2 


du'  dv'  . 
—  +  —  =  0 
dx  dt; 


where  the  second  equation  has  been  multiplied  by  v1/z.  Assuming  that  the  transformed 
field  variables  can  be  expanded  according  to 


u'(x,^;t)  =  u0(x,^;t)  +  vul(x,i^;t)  +  O(v2), 
v'(x,^-,t)  =  v0(x,i;-,t)  +  wl(x,^;t)  +  O(v2),  (145) 

p'  (x,  C f)  =  p0  (x,  C,  ;i t)  +  vp,  (x,  £ t)  +  0(v 2 ) , 


the  equations  for  the  components  are  found  by  collecting  terms  with  the  same  power  of  v 
as  a  factor.  The  zero-order  component  system  is  found  to  be 


0(1) 


«0 


du0 

+  «o-T2-  +  vo 
OX 


3“o 


+ 


dx  dC,1 


<K  |  _  o 

dx  dC 


(146) 


0(1) 

Bt  [w0  (x,  ; ;  t),  v0  (at,  £;t),p0  (x.  £ ;  r)] = Pt 

for  condition  at  v  =  C  =  0  (147) 

«o(*,?iOU=«((*.A). 

0(1) 

v0(*>(;OU=v1/Jvi(avi/jO. 

(148) 

p0(x,;-,t)  U=a(*>v‘/20, 

The  equations  for  the  first-order  system  are 


du< 


du 


du , 


du0  dpi  d\  d\ 


0(v) 


K'+M^+w'17+Voac+v'lc  +  &-ac2 

v,+„,^+v„^+^=ii 


dx2 


dx 


dt;  d£  dC 

5m!_+5vl  =  0 

dx  dt; 


(149) 
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0(v) 

5i[«I(x,^;0,v1(x,C;0,^1(^C;0]=0  for  condition  at  y  =  £  =0 

(150) 

«i |,=0=0, 

G(v) 

*** 

o 

II 

o 

(151) 

Uo=0, 

The  zero-order  system  is  the  Prandtl  boundary  layer  equations  with  the  addition  of  a  time 
varying  term  (these  equations  are  also  verified  by  the  Mathematica  program).  They  are 
completely  independent,  that  is,  they  only  involve  the  zero-order  field  variables.  However, 
the  first-order  system  is  coupled  to  the  zero-order  system.  Yet,  unlike  the  zero-order 
system,  the  first-order. system  is  linear.  Assuming  that  the  zero-order  system  can  be 
solved,  the  solution  can  be  inserted  into  Eq.  (149)  to  yield  a  linear  system  for  the  first-order 
field  variables. 

We  note  that  the  dynamics  equations  for  the  above  system  may  be  expressed  in  an 
altemate’form  after  some  manipulation.  The  alternate  forms  are  as  follows: 


8  (  2 )  S  /  \  dpo 

“•  +  &k  >+dS{U°V°)+-^  = 


d2u 


0(1) 


dx 
dPo 
Si 

Sup  |  dv0 
dx  8i 


Si 


=  0 


dr  \  s  (  \  8px  82u]  82u0 

\  ■  3  /  \  d  (  2  ^  dp\  d  v0 

0(v)  v,t_(„tv,)+-(v,  )+  — 


8u 

~8x  8i 


i  Sv, 

1  +— '-  =  0 


(152) 


(153) 


4.2. 1.3  Matching 

Once  the  solutions  to  the  two  separate  expansions  are  found,  they  must  be  “matched”  to 
one  another  to  deliver  the  full  composite  solution  across  the  entire  domain.  The  matching 
process  may  be  though  of  as  the  imposition  of  boundary  conditions  in  the  region  of 
common  validity.  For  the  zero-order  solutions,  the  process  is  fairly  simple.  The  higher 
order  expansions  require  more  effort. 

The  zero-order  outer  layer  solution  simply  is  the  Euler  equations.  In  the  zero-order 
approximation  the  Euler  solution  is  valid  everywhere  in  the  fluid  domain  except  for  an 
infinitesimally  thin  region  around  the  plate,  specifically,  a  thin  neighborhood  around  y=  0. 
The  inner  layer  expansion  is  valid  in  the  entire  region  of  the  transformed  coordinate 
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However,  since  C^ylv12  the  limit  as  v  approaches  zero  requires  that  C,  approach  oo  (for 
finite  y).  Thus,  our  zero-order  matching  conditions  are  given  by 


lim  H0  (*,  £ ;  t)  =  lim  U0  (x,  y;  t) = U0  (x,0;  t), 

£->oc  _V->0 

lim  v0  {x,  £ ;  t)  =  lim  v : l,2V0  (x,  y;  t) - v W2V0  (x,0;  t ),  (154) 

£  — >oo  V-»0 

lim  p0  (x,  Q\t)  =  \\mP(j  (x,  y;t)=P0  (x,0;  t). 

£  ->QC  y->0 

Matching  of  the  first-order  solution  is  more  involved.  The  two  solutions  are  required  to 
agree  to  first  order  over  an  open  interval  of  common  validity,  say  ys(y^,y2).  To  do  this 
typically  one  picks  an  intermediate  variable  of  the  form  yn=y/vn  and  transforms  the  two 
opposing  solutions  to  this  variable.  The  solutions  are  then  required  to  agree  to  first  order 
for  r)  in  (0,1). 

4.2.2  Example  -  The  Flat  Plate  Problem 

To  demonstrate  the  use  of  the  above  results  we  compute  the  steady-state  asymptotic 
solution  to  our  infinitesimally  thin  flat  plate.  Referring  to  Figure  10,  the  plate  is  placed 
edgewise  in  a  uniform  flow  at  the  position  y=0.  It  is  well  known  that  there  exists  a  semi- 
analytic  solution  for  the  zero-order  system,  the  Blasius  solution  [8].  Here,  however,  we 
add  the  first-order  corrections. 

4.2.2. 1  The  Zero-Order  Solution 

Due  to  the  presence  of  viscosity,  no-slip  boundary  conditions  are  imposed  upon  the  fluid 
at  the  plate  boundary.  At  the  edge  of  the  boundary  layer  (i.e.,  the  limit  as  C,  approaches  =o) 
the  fluid  must  assume  the  uniform  flow  conditions.  The  outer  layer  solution  is  simply  the 
solution  to  the  Euler  equations,  in  this  case  the  uniform  flow 

U0  =  ue(x,  y:  t)  =  =  const 

V0  =ve(x,y,t)  =  0  (155) 

P0=Pe(x,y;t)  =  Px=  const 

(As  mentioned,  the  subscript  e  indicates  the  solution  to  the  Euler  equations.)  The  zero- 
order  matching  conditions  require  that 

lim  uJx,C;t)  =  lim  ue(x,y,  t)  =  Ux, 

£  — >oo  y-¥  0 

lim  v0  (x,  £;  t)  =  lim  v  1/2ve(jc,  y;  t)  =  0 ,  (156) 

£  — >oo  y-*0 

lim  p0(x,<Z;t )  =  lim  pe(x,y;t)  =  Pxl . 

£  — >oo  y->0 

Along  with  these  matching  conditions,  we  also  have  the  following  boundary  condition  on 
our  asym  ptotic  system : 


u0(x,y\t)=Ua> 
v0(x,y;t)  =  0 


for  x  <  0 


u0(x,O;t)  =  O 
v0(*,0;f)  =  0 


for  x  >  0 


(157) 


In  the  region  x>0,  the  steady  state  solution  to  the  zero-order  boundary  layer  is  given  by  the 
Blasius  solution  [8],  which  is 
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«„(*,£)  =  tfx/'(» 7(x.O) 

v»(x,0  =  ~[^’(»7)-/(l)I 
PoM  =  P* 


rj(x,Q^ 


where  the  function  /fa)  solves  the  differential  system 


£f_+Lf£L 

drt !3  2  dr]1 


=  0 


/(0)  =  0 
/'(  0)  =  0 
lim/’(J?)  =  l 


(159) 


In  the  region  x<0  the  solution 
is  simply  the  uniform  flow 
conditions.  Thus,  the  fluid 
does  not  “feel”  the  plate  until 
it  reaches  the  region  x>0. 
Physically,  we  expect  a  wake 
extending  into  the  region  x<0 
around  the  front  edge  of  the 
plate,  the  point  (x,y)=(0,0). 
However,  our  analysis  does 
not  account  for  this  and  it  is 
well  known  that  the  Blasius 
solution  is  not  valid  in  a 
neighborhood  of  (0,0)  -  in 
fact,  it  is  singular. 


•  Figure  11:  The  Blasius  function  /  and  its  first 
derivative 


We  list  a  few  properties  of  the  function  /fa)  for  later  use.  It  can  be  shown  that  the 
boundary  conditions  on  /fa)  imply  [8] 


/"(0)  =  a  »  0.33206.  (160) 


We  also  have  the  useful  facts  that 


lim  (77  -f(ri))=  p  « 1.7208,  (161) 


and 


/"(*?)- 


(162) 


Finally,  the  first  few  terms  in  the  Taylor  expansion  are 


77*  a 2  r]5  lla  77 

/(77)  =  a-i - —  + - --  +  ... 

2!  2  5!  4  8! 


(163) 
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Using  the  initial  values  for  f{r\)  and  f’\r\) 
the  function  f(r\)  may  be  easily  computed 
numerically.  A  plot  of  f{r\)  and  f\r\)  is 
shown  in  Figure  11.  Note  that  /(r|) 
approaches  its  limiting  value  around  ri»5. 
The  streamlines  for  the  Blasius  solution 
for  LL- 1  are  shown  in  Figure  1 2.  These 
are  simply  the  level  sets  of  the  stream 
function  y=Uxf{r\)  and,  consequently,  the 
level  sets  of  r\(x,Q.  Due  to  the 
discontinuity  in  the  plate  at  x=0  we  can 
see  how  the  solution  is  invalid  in  the 
neighborhood  of  (x,y)=0.  Collecting  these 
results  we  find  that  the  zero-order  composite 


•  Figure  12:  Streamlines  for 
Blasius  solution 


solution  for  the  flat  plate  problem  is  given  by 


u(x,y)  =  Ux 

v(x,  y)  =  0  for  x  <  0 

p(x,y)  =  P* 


u(x,y)-UtBf' 


U  v 
2  x 

p(x,y)  =  P,a 


/ 

1  \ 

fu  } 

2 

y 

V*  j 

V 

) 

i  A 


k  w  / 


£/„ 


v  *  ) 


A 


fU.  ^ 


1 


\  wr  , 


for  x  >  0 


(164) 


4.2.2.2  The  First-Order  Solution 

Substituting  the  zero-order  outer  solution  into  the  first-order  outer  layer  equations  yields 
the  system  we  must  solve  for  Ui,  1/i  and  Pi.  Since  the  Euler  solution  is  simply  the  free¬ 
flow  conditions  Uo=Ux,  V0=0,  P0=P»,  we  have 


U, 


dU ,  dP, 


+  •— 4-  =  0 


dx  dx 
dx  dy 


=  0 


at/,  ^ 

dx  dy 


=  0 


(165) 


These  equations  can  be  combined  and  manipulated  to  yield  the  following  results: 


V2Vl(x,y)  =  0  (166) 

V2Pi(x,y)  =  0 


where  V2  denotes  the  Laplacian  operator.  The  solution  to  the  Laplacian  equations  can  be 
found  using  separation  of  variables.  After  imposing  the  boundary  conditions  at  x=0  and 
y-*»  we  find 
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V{  (x,  y)=\  Cj,  (i)  sin  foe  dk 

. . . Jl  . .  ,-;.r  . .  :  ^  (167) 

Px  (.V,  >’)  =  Cp (k)  e'A>'  sinkxdk 


where  Cp(k)  and  Cfk)  are  arbitrary  functions  of  the  separation  parameter  k.  Actually,  from 
inspection  we  see  that  Cfk)  is  the  Fourier  sine  transform  of  the  function  Pi(x,0),  likewise 
for  Cv(k)  and  V^x,0).  Thus,  these  functions  are  determined  from  the  matching  conditions 
between  the  two  layers. 


Now  consider  the  inner  layer.  Referring  to  the  first-order  inner  layer  equations  (the  O(v) 
equations),  the  second  equation  in  that  set  contains  expression  for  Uo,  v0  and  pi  only;  that 
is,  pi  is  the  only  unknown.  By  substituting  the  Blasius  solution  for  uq  and  Vq  into  this 
equation  we  get  the  following  relation: 


^L  =  if^>f7(/'(r7))2  +2UJ"(r1)-UJ(n)f\ri)]  (168) 

K  H  x  J 

*  % 

It  is  possible  to  integrate  this  expression  analytically  to  obtain  an  expression  for  Pi(x,Q. 
First  we  define  the  function  g(r|)  by 


g{f)  =  (169) 

0 

We  can  find  the  expression  for  g(r\)  by  integrating  the  differential  equation  for  /(r|).  Upon 
integrating  the  second  term  by  parts  we  find 

g(v)  =  2f’(T])  +  f(r])f'(T])-2a  (170) 

With  this  relation  one  may  determine  that  (again  using  integration  by  parts) 

fVj/'fo-)]2 da  =  [rig(T])l  -[ g{a)dc r 

Jo  0  (171) 

= 2nT07)  +  nf{ri)nri)  -  2/'(r?) 

Thus,  the  expression  forpi  is  given  by 

P\  (*,  £ )  =  7— IW"  (r?)  +  nf  (P)/'(P)  -  if  (n)f  ]+  ci  (*)  (172> 

4  x 

where  C)(x)  is  the  constant  of  integration  with  respect  to  0  To  find  Ci(x)  it  is  necessary  to 
apply  the  boundary  condition  at  x=0,  namely  pi  (0,0=0. 
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(173) 


lim  p]  (x,0  =  lim  777*'’’  +  lim  777  f(v)hf'  (v)  ~  fin)]  +  C,  (ux£ 2  /  rf ) 

*->o  2  Q  t)->oo4^z 

=  -  mb  C^/7 r) 

n-+*>  4  £ 

=  umlnLpm  +  q(u^2/r12) 

=  0 

This  implies  that  Ci(x)  is  of  the  form 

Ci(x)  =  -^f{juj^)+0(xn  (174) 

where  y>0.  With  such  an  expression  for  Ci(x)  the  pressure  is  zero  at  x=0,  except  at  the 
point  (x,0=(0,0)  where  the  Blasius  solution  is  known  to  be  singular. 

Next  we  attempt  to  match  the  inner  and  outer  solutions  for  the  pressure.  To  do  this  we 
introduce  the  variable 


(175) 

where  ye (0,1/2).  Converting  both  the  full  (i.e.,  up  to  order  0(v))  inner  and  outer 
expressions  to  this  variable,  they  should  agree  to  first  order  for  y  in  the  interval  (0,1/2)  as  v 
approaches  zero.  Letting 


1 


(176) 


we  get 


Px  (x,  y)  =  P„  +  v  7  ■ —  W'  (4) + 4f  (4)/'  (4)  -  (/(4))2  ]+  vCt  (*) 


4  x 

aP  +vLE^pm+vCl(X) 

4  x 


(177) 


\  X  / 


p2y+v\ 


4  x 


and 


P,  (x, y)  =  Paa  +vf  CP (k) e  kvr'’  sin  kx dk 

Jo  ^  (178) 

^  P_r  +v J  CP (k)  smkxdk-vl+ryj^Cp(k)ksmkxdk 

Because  of  the  CXy12*1)  term  in  the  inner  layer  solution,  it  is  impossible  to  match  the 
solutions  to  first  order.  This  fact  suggests  that  our  original  assumption  that  the  solutions 
could  be  expanded  in  powers  of  v  was  probably  wrong.  A  more  likely  expansion  would 
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then  seem  to  be  in  powers  of  v1/2.  However,  due  to  time  constraints  in  Phase  I  we  have 
■  not  been  able  to  further  pursue  this  avenue  of  investigation. 

4.2.3  Wall  Curvature 

Here  we  briefly  discuss  the  situation  were  the  wall  is  curved,  rather  than  a  simple  flat  plate. 
Assume  that  we  have  a  coordinate  system  as  in  Figure  9,  where  the  tangential  coordinate 
is  given  by  x.  We  denote  the  curvature  of  the  wall  by  k(x),  where  k(x)  is  assumed  to  be 
continuous  and  have  a  first  derivative.  Using  a  boundary  layer  analysis  where  Q  is  the 
stretched  normal  coordinate  as  before,  it  can  be  shown  that  the  curvature  will  affect  the 
pressure  distribution  according  to  the  following  [8]: 

^=v1,2k(x)u'2  (179) 

In  our  present  expansion  of  p’  this  equation  would  probably  have  to  be  added  into  the  0(1) 
system  and  would  affect  the  boundary  layer  directly.  However,  if  we  tried  an  expansion  in 
powers  ®f  v1/2  the  curvature  would  enter  into  the  0(v1/2)  system.  From  the  previous 
analysis  we  expect  that  system  to  be  linear  in  the  field  variables.  Thus,  it  may  be  possible 
to  incorporate  wall  curvature  and  achieve  some  sort  of  analytic  predictions.  This  condition 
would  be  extremely  useful  if  one  were  to  employ  the  curvature  as  the  control,  for  example, 
using  MEMS  type  devices  mounted  on  the  body  surface.  As  with  the  case  above,  our 
funding  for  Phase  I  ran  out  and  we  were  not  able  to  pursue  this  idea  as  of  yet. 

4.2.4  Multiple  Time  Scales 

The  time  derivative  of  Vq  appears  in  the  first-order  equations  and  not  the  zero-order 
equations.  However,  v0  and  u0  are  not  independent,  they  are  related  by  the  divergence 
conditions.  Thus,  the  time  behavior  of  v0  will  be  dictated  by  the  time  behavior  of  Uo  through 
that  condition.  Yet  it  is  interesting  to  investigate  the  possibility  of  two  time  scales,  a  slow 
time  scale  f,  and  a  fast  time  scale  t2.  Such  a  condition  would  be  necessary  to  introduce 
the  time  derivative  of  vb  into  the  zero-order  system,  which  are  the  boundary  layer 
equations.  The  form  of  Eqs.  (144)  suggests  the  following  definitions 

tx=t,  and  r,  =  — ,  (180) 

"  v 


so  that 


JL-jL  iJL 

dt  dti  v  dt2  ’ 


(181) 


and  the  field  variables  are  now  functions  of  x,y,fi,  and  t2  (e.g.,  i/(x,y,U,t2)  )•  Applying  this 
transform  to  Eqs.  (144)  yield  the  system 
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(182) 


du'  1  du'  ,du'  ,du'  dp' _  d2u'  d2u' 

dr,  v  dt2  dx  dQ  dx  dx2  dC,2 


dv'  dv' 

V  —  +  —  +  v 
dt{  dt2 


,dv'  ,dv') 
u —  +  v' — 
dx  dC,  j 


d2v' 


d2v' 


dp'  2 

+  ^—=v  — r  +  v  , 
dt;  dx2  dQ2 


du'  dv' 
-fo+~d<Z 


=  0 


Inserting  the  expansion  of  Eqs.  (145)  into  the  above  yields  the  following  systems  of 
asymptotic  expressions 


0(1 /v) 


(183) 


0(1) 


duL+du1 
dt2  dtl 


du0  dun  dp 0  <32u0 

+  M„ - -  +  V0 - -  +  -^  = - f 

0  dx  dC  dx  dQ- 

|  <?Po  _  0 
dt2  dt; 


(184) 


duL+dv1 

dx 


=  0 


Qy) 


du{ 

dt 


■  +  Ur, 


dux 

dx 


-  +  w. 


duo 

dx 


■  +  vn 


du* 


-  +  v. 


ou . 


di 


d\  d\ 

dt;1  dx2 


dv,  dv0  3v0 

dt2  dtx  dx 


+  vo 


dv0  |  qgi  _  d\ 
dt;  dt;  dt;2 


(185) 


dux  dv, 
dx  dC ' 


=  0 


There  also  exist  boundary  conditions  associated  with  the  above  system.  This  situation  is 
more  complicated  than  the  single  time  case  since  now  we  have  coupling  between  Uo  and 
u,  in  the  zero-order  (  0(1) )  equations.  However,  it  is  possible  to  decouple  the  system  by 
employing  the  Fredholm  alternative  theorem  [7]. 

Note  that  the  initial  conditions  in  Eq.  (184)  require  that  Lh(x,y,t i,fe)=0  when  £=fi=fe= 0.  We 
now  impose  the  additional  constraint  that  ui=0  at  f2=77v,  where  we  shall  take  t=T as  the 
final  time.  Thus, 

II, (x,  y\tw)  eX  =  {he  C'([0,77v))|  h(  0)  =  0,  h(T/v)  =  0  J.  (186) 

For  clarity  let  L  denote  the  operator  d/dt2  and  define  the  function  space  Y according  to 
Y  =  \he  C°([0,77v))}  (187) 

so  that 

L:X->Y 

h»dh/dt2  °88) 
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Both  Xand  / can  be  given  a  Hilbert  space  structure  by  introducing  the  inner  product 


(g,^)sfg( tiyWi)dt2  (189) 

0 


Then  the  adjoint  operator  L*  is  easily  found  to  be 


L':X->Y 
h  1 — )  —  dhj  dt2 


(190) 


using  a  simple 
integration  by  parts. 

From  the  lack  of 
boundary  conditions  in 
the  space  Y  we  see  that 
L*  has  a  nontrivial  null 
space."  To  insure  that 
the  first  of  Eqs.  (184)  has 
a  solution,  the  Fredholm 
alternative  requires  that 
the  terms  involving  the 
zero-order  components 
must  be  in  the  range  of 
L=8/5t2.  Recall  from 
Hilbert  space  theory  that 
the  range  of  the  operator 
L  is  always  orthogonal  to 
the  null  space  of  the  adjoint  L*.  That  is, 


•  Figure  13:  Schematic  of  Fredholm  alternative 


Y  z>  Range(L)-LNull(L  )  c:  Y.  (191) 

This  condition  is  represented  graphically  in  Figure  1 3  where  the  only  intersection  between 
Range(L)  and  Null{L*)  is  the  zero  vector  (likewise  for  RangeiL*)  and  Null(L)).  Since  the 
null  space  of  L*=dl8t2  is  the  set  of  constant  functions  on  [0,77v),  the  orthogonality  relation 
requires  that 


Tj  v 


Sup 

dt\ 


+  un 


<K 

dx 


'  +  Vn 


du0  fop 
dx 


d2u0 

dt2 


\dt2=  0 


J 


(192) 


for  all  x,  y,  and  f|.  Practically,  this  condition  is  met  only  when  the  integrand  is  identically 
zero  (another  alternative  is  that  the  integrand  is  periodic  in  t2,  an  unlikely  event). 
Therefore,  5^2=0  so  that  Ui  is  also  independent  of  the  fast  time  scale.  Applying  a 
similar  analysis  to  the  second  of  Eqs.  (184)  implies  that  cVi/<3?2=0.  Collecting  these  results 
we  have 
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Temporal  boundary  conditions: 


Uq  (x,  y;0,0)  =  U,  (x,  y)  v0(x,  y;0,0)  =  v, (x,  y)  p0  (x,  y,0,0)  =  p,  (x,  y) 
u,(x,y;0,0)  =  0  v,(x,y;0,0)  =  0  p,(x,>';0)0)  =  0 

u,(x,y;t,T/v)  =  0  vl(x,y;t,T/v)  =  0 

Matching  conditions: 


(193) 


lim  u0  (x,  £ ;  t)  =  lim  ue(x,y;t), 


v->0 


limv0(x,£;O  =  limv  ve(x,y;t ), 

£  — >oc  V->0 

lim  p0(x,C;t)  =  lhnpe(x,y,t). 

£  -»oo  y— >0 

-  f 


(194) 


Independence  of  fast  time  scale: 


dun  „  3ul=0  gv) 


Sr2  gr2  gr2 


=  0 


(195) 


Dynamics  equations  0(1): 


5wo 

gr. 


•  +  «„ 


g“o 

dx 


■  +  Vn 


du0 

,  Spo  . 

d2u{ 

dx 

gv0 

,dPo  . 

=  0 

dt2 

dt 

g«o 

-O 

dx 

—  \J 

(196) 


Dynamics  Equations  0(v): 


g«, 

gr, 


•  +  m. 


g«, 

5x 


+  M, 


gw 

dx 

5v o 
gr. 


-  + v„ 


+  M„ 


ac 

gy0 

gx 


■  +  V, 


g«n 


g/>,  _  g2«. 


,  d2»o 

g£  g*  g<?2  3x2 


-  +  Vft 


gv„ 


g/>i  .  g2Vp 


g£  gf  g£2 

5«,+5vl  =  0 

gx  g£ 


(197) 


These  equations  could  provide  useful  when  attempting  to  design  a  controller  where  the 
response  time  is  comparable  to  the  flow  rate  reaction.  This  possibility  was  mentioned 
when  discussing  the  control  of  the  discrete  case. 
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4.3  Conclusion 

As  we  have  mentioned  we  are  still  investigating  the  control  aspects  of  this  project  and  our 
approach  so  fir  has  been  rather  pragmatic.  It  is  hoped  that  the  foray  into  the  theoretical 
structure  of  the  Navier-Stokes  equations  can  provide  some  practical  insight  into  fluid-flow 
controller  design.  This  is  especially  true  with  respect  to  the  last  sections  where  we  hope  to 
determine  the  dynamics  of  active  boundary  layer  stabilization. 
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APPENDIX  A:  Source  Listing 


The  Mathematica  source  file  used  in  the  boundary  layer  analysis  is  listed  here.  The 
source  symbolically  computes  the  two-dimensional  incompressible  Navier-Stokes 
equations  for  any  coordinate  system.  The  general  coordinates  are  (u1,u2)  and  their 
mapping  into  cartesian  space  is  supplied  by  the  functions  x1(u1,u2)  and  x2(u1,u2).  The 
principles  of  differential  geometry  are  then  applied  to  generation  the  form  of  the  Navier- 
Stokes  equations  in  the  general  coordinates  (u1,u2).  The  current  listing  is  for  stretched 
polar  coordinates  using  a  stretching  parameter  a.  Specifically,  the  radial  coordinate  is 
stretched  around  the  value  u1=1  by  a  factor  1/a.  At  the  end  of  the  file  the  stretching 
parameter  a  is  taken  to  be  v2  and  limits  are  taken  for  the  regime  a->0.  The  boundary 
layer  equations  are  then  recovered. 
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GeneralFluidFlowWithStreching.nb 


1 


Fluid  Flow  in  Polar  Coordinates 

■  Preliminaries 

ClearAll[] 

■  Load  Packages 

<<  "G:\Wolfram  Resear ch\Mathematica3 . 0\AddCtas\StandardPackages\Graphics\PlotField.m" 
Names  [  "Graphics '  PlotField'  *  "  ] 

{} 

■  Our  favorite  standard  basis 

el  =  {1,  0} 

{1.  0} 

e2  =  {0,  1} 

{0,  1} 


■  Stretched  Polar  Coordinates 

We  expand  the  traditional  polar  coordinate  r  about  r=l  by  the  factor  a. 
Our  new  coordinate  being  denoted  by  £  the  result  is  (r-  l)/a.  OR 
r  =  1  +  ar£ 

The  Cartesian  map  is  then  given  by 

x(£0)  -x\  (£0)ei+r2(£  0)e2 

=  (l+a£)cos#  ej  +  (l+or£)sini9  e2 

■  The  general  coordinates  -  (ul,u2) 

Clear  [ul,  u2] 
u  =  {ul,  u2} 

{ul,  u2} 
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■  Set  up  the  Cartesian  map  x[ul,u2]  =  {xl[ul,u2],  x2[ul,u2]} 

xl[ul_,  u2_J  =  (1  +  aul)  Cos[u2] 

(1  +  ul  a)  Cos [u2] 

x2  [ul _ /  u2_]  =  (1  +  aul)  Sin[u2] 

(1  +  ul  a)  Sin[u2] 

x [ul _ /  u2_]  =xl[ul,  u2]  el  +  x2[ul,  u2]  e2 

{ (1  +  ula)  Cos[u2],  (1  +  ula)  Sin[u2]} 


■  Set  up  the  tangent  plane  basis  vectors:  T(Ui)U2)  R  -  span{gl,g2} 

gl[ul_,  u2_J  x[ul,  u2] 

{a  Cos [u2] ,  a Sin[u2] } 

g2  [ul_,  u2_]  =  3u2  x[ul,  u2] 

{-(1  +  ula)  Sin[u2],  (1  +  ula)  Cos[u2]} 


■  Now  compute  the  Jacobian  and  the  dual  basis  for  the  Cotangent  bundle  (grl  and  gr2) 

G[ul_,  u2_]  =  Transpose [{gl[ul,  u2] ,  g2[ul,  u2]}] 

{ {a  Cos  [u2],  -(1  +  ula)  Sin[u2]},  {aSin[u2](  (1  +  ula)  Cos[u2]}} 


Ginv[ul_,  u2_]  =  Simplify [  Inverse  [G[ul,  u2]  ]  ] 


r r  COS [U2] 
^  a 


Sin[u2]  ■( 

5  J' 


r  Sin[u2] 
i  1  +  ula  ' 


Cos [u2]  n 
1  +  ul  a  J  J 


MatrixForm[G[ul,  u2]  ] 

/aCos[u2]  -(1  +  ula)  Sin[u2]  \ 
\  a  Sin  [u2  ]  (1  +  ul  a)  Cos  [u2]  / 


MatrixForra[Ginv[ul,  u2]  ] 

r  Cos  fu2]  Sinfu2]  ' 

a  a 

_  Sinru21  Cos fu21 

t  1+ula  1+ula 


J[ul_,  u2 U  =  Simplify [Det [ G [ul ,  u2]  ]  ] 

a  (1  +  ul  a) 


grl [ul_,  u2 _]  =  Ginv[ul,  u2]  [  [1]  1 


j-  Cos  [u2] 


Sin[u2] 


} 


a 


a 
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'  _  '  gr2 [ul_,  u2_]  =  Ginv[ul,  u2]  [  [2]  ] 

r  Sin[u2]  Cos[u2]  ■. 
t-  1  +  ul  a  *  1  +  ul  a  ■* 

■  The  metric  induced  from  E2  -  denoted  g 

,  „•  , r/gl[ul,  u2]  .gl[ul,  u2]  gl[ul,  u2]  .g2[ul,  u2]  v, 

g[ul_ ,  u2_J  =Sin®llfyL(g2[ui/ u2]  .gl[ul,  u2]  g2[ul,  u2]  .  g2[ul,  u2]  /  J 

{{a2,  0},  {0,  (1  +  ula)2}} 

MatrixForm[g[ul,  u2]  ] 

I  a2  0 

i  0  (1  +  ula)2 


■  The  contravariant  metric  (metric  for  covectors  in  basis  grl,gr2) 
ginv[ul_,  u2_J  =  Inverse[g[ul,  u2]  ] 

General: : spelll  : 

Possible  spelling  error:  new  symbol  name  "ginv"  is  similar  to  existing  symbol  "Ginv” . 

rr  (1+ula)2  Qi  rQ  _ a? _ ii 

>•  a2  +  2  ul  a3  +  ul2  a4  '  1 '  '  a2  +  2  ul  a3  +  ul2  a4 

MatrixForm[ginv[ul,  u2]] 

f  (1+ula)2  q  1 

a2 +2  ul  ct^+ul2  a4 

n  a2 

k  u  a2 +2  ul  a3 +ul2  a4  j 


■  Christoffel  Symbols 

We  setup  the  symbols  so  that  r|[ul,  u2]  =  T[ul,  u2][[i,  j,  k]] 

rl  [ul _ ,  u2_]  =  Transpose [ 

Simplify[ 

Ginv(ul,  u2]  .D[G[ul,  u2] ,  ul] ] ] 

(t°'  °>' t0'  rviinr}} 


r2[ul_,  u2_]  =  Transpose [ 

Simplify! 

Ginv[ul,  u2]  .  D[G[ul,  u2]„  u2]]] 


{{°< 


a 


1  +  ul 


1  +  ul  a 
a 


^  o}} 


r [ul_,  u2_J  =  Join[{rl[ul,  u2]},  (r2[ul,  u2]}] 

{{to.  0).  {0.  TTula'))'  {{”•  TT^nr).  »))) 
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u  Check  them 

r[ul,  u2]  [  [1/  2,  2]] 

a 

1  +  ul  a 

r[ul,  u2]  [  [2/  1,  2]] 

a 

1  +  ul  a 

r[ul,  u2]  [  [2,  2,  1]] 

_  1  +  ula 
a 

■  View  the  Constant  Coordinate  Lines  in  R2 

Well  do  this  for  a~Q2  and  (ul,n2)e[0,l]x[0,^/2] 

a  =  0 .2 

0.2 

tblLinesl[ul_]  =  Table[x[ul,  i0.1Pi/2],  {i,  0,  10}] 
ParametricPlot [Evaluate [ tblLinesl [ul]  ]  ,  {ul,  0,  1}] 


-  Graphics  - 

tblLines2[u2_]  *  Table[x[i  0.1,  u2Pi/2],  {i,  0,  10}] 
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1  ParametricPlot [Evaluate [ tblLines2 [u2]  ] ,  {u2,  0,  1}] 


-  Graphics  - 

tblLines  [x_]  =  Join[tblLinesl[x] ,  tblLines2[x]  ] ; 
ParametricPlot  [Evaluate  [  tblLines  [x]  ] ,  (x,  0,  1}] 


-  Graphics  - 


■  View  the  Polar  Tangent  Bundle 

tblTangl  =  Table[{x[0.1  i,  Pi/20j],  gl[0.1i,  Pi/20  j]}, 
{i»  0,  10}/  { j /  0/  10}]; 

tblTang2  =  Tablet {x[0.1  i,  Pi/20  j] ,  g2[0.1i,  Pi/20  j]}, 
{i,  0,  10},  {j,  0,  10}]; 
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■  tblTang  =  Flatten  [Join  [tblTangl,  tblTang2] ,  1]; 

Lis tPlotVectorFieldf  tblTang,  VectorHeads ->  True,  ScaleFactor  ->  0 . 1] 


-  Graphics  - 

■  Now  we  have  to  clear  off  a  for  symbolic  computation 
g2[ul,  u2] 

{  -  (1  +  0 .2  ul)  Sin[u2]  ,  (1  +  0.2  ul)  Cos  [u2]  } 

Clear[a] 
g2[ul,  u2] 

{-  (1  +  ul  a)  Sin[u2],  (1  +  ula)  Cos[u2]} 


■  Conversion  Between  Covectors,  Vectors,  and  Physical  Components 


■  Covectors  to  Vectors 

Sub2Sup[v :  = 

Simplify! 

{Sum[ginv[ul,  u2]  [  [1,  i]  ]  v[  [i]  ] ,  {i,  2>] , 
Sum[ginv[ul,  u2]  [[2,  i]]  v[[i]] ,  (i,  2}]} 

] 
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■  Vectors  to  Covectors 


Sup2Sub  [v :  = 

Simplify! 

{Sum[g[ul,  u2]  [  [1,  i]  ]  v[  [i]  ] ,  {i,  2}] , 
Sum[g[ul,  u2]  [[2,  i]]  v[[i]],  (i,  2}]} 

] 


■  Vectors  to  Physical  Components 

Sup2Phys[v  :  = 

Simplify[ 

{v[  [1]  ]  /  V ginv[ul,  u2]  [[1,  1]]  , 
v[  [2]  ]  /  VginvCul,  u2]  [[2,  2]]  } 

] 


■  Covectors  to  Physical  Components 

Sub2Phys[v  :  = 

Simplify [ 

{v[[l]]/Vg[ul,  u2]  [[1,  1]]  , 
v[[2]]/Vg[ul,  u2]  [[2,  2]]  } 

] 

General :: spelll  : 

Possible  spelling  error:  new  symbol  name  "Sub2Phys"  is  similar  to  existing  symbol  "Sup2Phys". 


■  The  Covariant  Derivative 


■  The  Covariant  Derivative  of  (Contravariant)  Vectors 

dCovSup[i,  v]  =  (VjV1,  VjV2} 

dCovSup[i_,  v  s  _}]  s  = 

{D[v[  [1]  ] ,  u[  [i]  ]  ]  /  D[v[  [2]  ] ,  u[  [i]  ]  ] } 

+  {Sum[r  [ul,  u2]  [  [i,  j,  1]]  {j,  2}], 

Sum[r[ul,  u2]  [[i,  j,  2]]  v[[j]],  {j,  2}]} 

■  The  Covariant  Derivative  of  Covariant  Vectors  (Covectors) 


dCovSub[i,  v]  =  {Vj vi,  VjV2} 
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dCovSub[i_,  v  :  {_,  _}]  :  = 

{D[v[  [1]  ] ,  u[  [i]  ]  ]  /  D[v[[2]],  u[fi]]]} 

-  {Sum[r[ul,  u2][[i,  1,  k]]  v[[k]],  {k,  2}], 

Sum[r  [ul,  u2]  [  [i,  2,  k]  ]  v[  [k]  ] ,  {k,  2}] } 

General: sspelll  : 

Possible  spelling  error:  new  symbol  name  "dCovSub"  is  similar  to  existing  symbol  "dCovSup" . 


■  Define  Various  Calculus  Operators 

■  These  operators  take  scalars  to  covectors 

Grad[f_]  :=  {D[f ,  u[  [1]  ]  ] ,  D[f ,  u[  [2]  ]  ] } 

■  These  operators  take  scalars  to  vectors 

Curl2Vector  [  f_J  :  = 

Simplify! 

(1/  J[ul,  u2] )  (D[f,  u2]  ,  -D[f ,  ul]} 

1 

■  These  operators  take  covectors  to  scalars 

Curl2  Scalar  [v  :  = 

Simplify! 

(1/J[ul,  u2]) 

(dCovSub[l,  v]  .  e2  -  dCovSub[2,  v]  .  el) 

1 

■  These  operators  take  vectors  to  scalars 

Diverg[v  :  {_,  _}]  :=  dCovSup[l/  v]  .  el 
+  dCovSup[2,  v]  .  e2 

■  These  operators  take  scalars  to  scalars 

T,ap] a<-i angpal ar[fj  :=  Sum[  ginv[ul,  u2][[i]]  ,dCovSub[i,  Grad[f]],  {if  2}  ] 

■  These  operators  take  vectors  to  vectors 

GradSup[f_]  :=  Module!  {df},  df  =  Grad[f]  ? 

Simplify! 

{Sum[  ginvful,  u2]  [  [i,  1]  ]  df  [  [i]] ,  (i,  2}] , 

Sum[ginv[ul,  u2]  [  [i,  2]  ]  df  [  [i]  ] ,  {i,  2}] } 

1 
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’Mvec[v  :  {_,  _}]  :=  v[  [1]  ]  dCovSupfl,  v] 
+  v[  [2]  ]  dCovSup[2,  v] 

LaplacianVector[v  :  _}]  :  = 

Simplify[ 

{Sum[ginv[ul,  u2]  [[j,  k]] 

(dCovSup[j, 
dCovSup[k,  v] 

]  •  el) ,  { 3 /  2} ,  {k,  2}], 
Sum[ginv[ul,  u2]  [[j,  k]  ] 

(dCovSup[j, 
dCovSup[k,  v] 

]  .e2),  {j,  2},  {k,  2}] 

} 

] 


■  Compute  the  Vector  Calculus  Operations  in  Polar  Coordinates 

Grad[f  [ul,  u2]  ] 

{f(1-0)  [ul,  u2]  ,  f,0-X)  [ul,  u2]  } 


GradSup[f [ul,  u2]] 

f  fi1’0*  [ul,  u2]  f10-11  [ul,  u2]  | 

L  a2  '  (1  +  ula)2  * 


Diverg[{vl[ul,  u2] ,  v2[ul,  u2]}] 

avl[ul,  u2]_  +  v2(o.i>  [ul,  u2]  +vl11,0’  [ul,  u2] 
1  +  ul  a 

Curl2 Scalar [ (vl [ul,  u2] ,  v2[ul,  u2]>] 

-vl10'11  [ul,  u2]  +  v2(1,01  [ul,  u2] 
a  +  ul  a2 


Curl2Vector[f  [ul,  u2]] 


{ 


f10-11  [ul,  u2] 
ot  +  ul  a2 


fd.Q)  [ul,  u2]  , 
a  (1  +  ul  a)  > 


Advec[{vl[ul,  u2] ,  v2[ul,  u2]>] 


{v2  [ul,  u2] 
v2 [ul,  u2] 


/  (1  +  ula)  v2[ul,  u2]  +vl(0,i) 

\  a 

lavUul,  u2]_+v2io.u  [ul/  u2]) 
\  1  +  ula  l 


[ul,  u2]  |  +vl[ul,  u2]  vla,0)  [ul,  u2] , 

+ vl[ul,  u2]  (£^i^2]_+v2a.o)[ul(  u2])} 


LaplacianScalar[f  [ul,  u2]] 


a2  (f  <°-2>  [ul,  u2]  +  (1*ulal£l*'°’  (ui,u2L)  (1  +  ula)2  f(2,0>  [ul,  u2] 
a3  +  2  ul  a3  +  ul2  a4  a2  +  2  ul  a3  +  ul2  a4 
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'  LaplacianVector[{vl[ul,  u2] ,  v2[ul,  u2]}] 

f - i - -  (-a2  vl[ul,  u2]  -2a  (1  +  ul a)  V2'0,11  [ul,  u2]  + 

1  a2  (1+ula)2 

a2  vl<0,2)  [ul,  u2]  +  vl<2,0)  [ul,  u2]  +2  ula  vl<2,0)  [ul,  u2]  +  ul2  a2  vl<2,0)  [ul,  u2] ) , 

- - - r  (-a2  (1  +  ul  a)  v2  [ul,  u2]  +  2  a3  vl(0,1)  [ul,  u2]  + 

a2  (1  +  ula)3 

(1+ula)  (a2 v2(0-21  [ul,  u2]  +  2av2a'0’  (ul,  u2]  +  2ula2v2a,0)  [ul,  u2]  + 
v2(2,0)  [ul,  u2]  +  2  ula  v2<2,0)  [ul,  u2]  +  ul2  a2  v2<2,0)  [ul,  u2] ) )  } 


Finally:  Fluid  Flow 

Flow  around  the. cylinder 

■  First  Find  Boundary  Conditions  for  Uniform  Flow 


■  Uniform  flow  in  x-direction  in  Cartesian 

VinfCart  =  {1,  0} 

{1,  0} 


■  Convert  to  Polar 

Vinf [ul_,  u2_]  =  Ginv[ul,  u2]  .VinfCart 

C  Cos  [u2]  Sin[u2]  -> 
i  a  '  ~  1+ula  1 


■  Find  Physical  Components 


Vinf Phys  [ul_,  u2_]  =  {Vinf[ul,  u2]  [  [1]  ] /Vginv[ul,  u2]  [[1,  1]  ]  , 

Vinf  [ul,  u2]  [  [2]  ]  /  V  ginv[ul,  u2]  [[2,  2]]  } 


Cos [u2] 


Sin[u2] 


(1+Ulg)2 

t2  ul  nr3  inl-2  0,4 


( 1  x  nl  nr ^ 


r 
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■  Take  limit  as  ul  approaches  0 


Limit  [VinfPhysful,  u2] ,  ul  ->  0] 


a  Cos [u2] ,  -Sin[u2] 


} 


Analytic  solution  to  inviscid  flow  around  a  unit  cylinder  (physical  components) 

1  \  1  1  Sin[u2] } 


Vcyl[ul_,  u2_]  =  {1 - -  Cos[u2],  -  1+  ■ 

l(  (1  +  a  ul)2  J  [  (1  +  a  ul) 

{ (l - - - j-]  Cos  [u2] ,  (-1-  — — \ — r —)  Sin[u2] } 

M.  (1  +  ul  a) 2  /  V  (1  +ul  a)  / 


-  » 


Pcyl[ul_,  u2_J  =  C  -  — 


(1  +  a  ul) 4  -  2  (1  +  a  ul) 2  Cos  [2  u2]  +1 
(1  +  a  ul)4 


General: :spelll  : 

Possible  spelling  error:  new  symbol  name  "Pcyl"  is  similar  to  existing  symbol  HVcyl  . 

1+  (1  +ula)4  -  2  (1  ul  a) 2  Cos  [2  u2] 

2  (1  +ul  a)4 


■  Take  limit  as  ul  approaches  0 

Limit  [Vcyl[ul,  u2] ,  a  ->  0] 

{0 ,  -2  Sin[u2] } 

Simplify  [Limit  [Pcyl  [ul,  u2] ,  a  ->  0]  ] 

-1  +  C  +  Cos  [2  u2] 


■  Define  the  Navier  Stokes  Operator  in  General  Coordinates 

Navi er Stokes  [v  P_/  vj  :  = 

Advec[v] 

-  y  LaplacianVector[v] 

+  Sub 2 Sup [  Grad[p]  ] 
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Do  the  Stretching  and  Take  Limits 


Define  our  starting  equations  with  parameters 

nsEql [ul_,  u2_,  v_J  =  NavierStokes[{vl[ul,  u2] ,  v2[ul,  u2]},  p[ul,  u2] ,  v] 


[v2  [ul,  u2]  |- 


(ltul  a)  v2[ul,  u2] 


+  vl10,11  [ul,  u2]  j 


,(1.0) 


[ul,  u2] 


+  vl[ul,  u2]  vla'0)  [ul,  u2]  - 

(v  (-a2  vl  [ul,  u2]  -  2a  (1  +  ula)  v2(0'1)  [ul,  u2]  +  a2  vl<0,2)  [ul,  u2]  + 


a2  (1  +  ul  a) 2 

vl'2'0>  [ul,  u2]  +  2  ul  a  vl(2'0)  [ul,  u2]  +  ul2  a2  vl(2'0)  [ul,  u2] ) ) , 
p10-11  [ul,  u2]  + 

(1+ula)2 

v2  [ul,  u2]  *v2,0,1>  [ul>  u2])  H-vl[ul,  u2]  +  v2a.0)  [ul>^  u2] )  . 

- - - r  (v  (-a2  (1  +ul  a)  v2  [ul,  u2]  +  2  a3  vl10'1’  [ul,  u2]  + 

a2  (1+ula)3 

(1+ula)  (a2  v2,0'2)  [ul,  u2]  +  2av2(1'0)  [ul,  u2]  +  2  ul  a2  v2(1'0)  [ul,  u2]  + 
v2<2’01  [ul,  u2]  +  2  ula  v2(2,0)  [ul,  u2]  +  ul2  a2  v2,2i0)  [ul,  u2] ) ) )  } 

nsEql  [ul_,  u2_,  v_J  =  Sup2Phys  [nsEql [ul,  u2,  v]  ] 

{^_  (v2[ul,  u2]  (-(lH-ulg)  ^^L+vi,o-»[ul,  u2] 


>U,0) 


[ul,  u2] 


a2 

1 


+  vl  [ul ,  u2]  vl(1'°»  [ul,  u2]  - 


—  (v  (-a2  vl  [ul,  u2]  -  2  a  (1  +  ula)  v2(0'1)  [ul,  u2]  +a2.vl(0>2)  [ul,  u2]  + 


a2  (1  +  ul  a) 2 

vi (2. °)  [ul(  U2]  +2  ul  a  vl(2,0)  [ul,  u2]  +ul2  a2  vl<2,0)  [ul,  u2] )  )| , 
1  (  pto-1’  [ul,  u2] 


v  a*uii 


(1  +  ula)2 
r  a  vl  [ul,  u2] 


v2  [ul,  u2]  (  aVi[+Uul^~  +v2,0,1)  tul'  u2l)  +vltul/  u2] 


a  v2 [ul,  u2] 
1  +  ul  a 


+  v2a,0)  [ul,  u2]  | 


_  (v  (-a2  (1  +  ula)  v2  [ul,  u2]  +  2  a3  vl(0'X)  [ul,  u2]  + 


a2  (1  +  ul  a) ' 

(1  +  ula)  (a2v210'2)  [ul,  u2]  +  2  av2ll'0)  [ul,  u2]  +  2  ul  a2  v2a‘0)  [ul,  u2]  + 
v2(2-0)  [ul,  u2]  +  2  ul  a  v2<2,0)  [ul,  u2]  +ul2  a2  v2(2,0)  [ul,  u2] ) ))  j } 
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f 

Now  multiply  radial  component  by  a 

nsEq2[ul_,  u2_,  v_J  =  {{a,  0},  {0,  1)}  .nsEql[ul,  u2,  v] 


a  |v2  [ul,  u2]  (- 
a'0)  [ul,  u2 


(1  +  ul  a)  v2  [ul,  u2] 


a 


+  vl(0ll)  [ul,  u2] )  + 


+  vl[ul,  u2]  vl(1'0>  [ul,  u2]  - 
—  (v  (-a2  vl  [ul,  u2]  -2  a  (1  +  ula)  v2(0,1)  [ul,  u2]  +  a2vl(0,2>  [ul,  u2]  + 


a2  (1  +  ul  a) 

vl(2-°>  [ul,  u2]  +  2  ul  a  vl<2,0)  [ul,  u2]  +  ul2  a2  vl<2'0)  [ul,  u2]) )  j| , 


^(0,1) 


[ul,  u2] 


j  nz  i  (i + ui a) 2 

V  (l+ul  a) 1 

v2[ul,  u2)  +  v2 10-11  [ul,  u2l)  ,vl[ul,  u2]  , v2a."  [ul,  a2) )  _ 

—  (v  (-a2  (1  +  ula)  v2[ul,  u2]  +  2  a3  vl'0,1’  [ul,  u2]  +  ' 


a2  (1  +  ul  a)  • 

(1  +  ul  a)  (a2  v2(0'2)  [ul,  u2]  +  2  av2(1,0)  [ul,  u2]  +  2  ul  a2  v2a-0)  [ul,  u2]  + 
v2,2,0)  [ul,  u2]  +  2  ul  a  v2(2,01  [ul,  u2]  +  ul2  a2  v2(2,0>  [ul,  u2] ) ) )  |  } 


Now  assume  that  the  stretching  parameter  or=Vv,  or  v  =  a2 


nsEq3[ul_,  u2_]  =  nsEq2[ul,  u2,  aA2] 

■  ul  a)  v2  [ul,  u2] 


{  -j==-  |a  |v2  [ul,  u2]  (-  —  — 

V  a2” 


h  vl10'11  [ul,  u2]  | 


V (1+ula) 


pa'°’  [U.1,  u2]-  +  vl  [ul,  u2]  vla'0)  [ul,  u2]  - 
a2 

_ l -  (-a2  vl  [ul ,  u21  -  2a  (1+ula)  V2'0'1’  [ul,  u2]  +  a2  vl(0'2)  [ul,  u2]  + 

(1+ula)2 

vl(2’0)  [ul,  u2]  +  2  ul  a  vl(2'01  [ul,  u2]  +  ul2  a2  vl<2,0)  [ul,  u2])||, 

p10-11  [ul,  u2]  + 

(1+ula)2 


v2  [ul,  u2]  (aVi[,Uu1iaU2j-+v2(°,1)  [ul,  u2] )  +  vl  [ul,  u2]  +  v2<^>  [ul,  u2; 

(-a2  (1  +  ula)  v2  [ul,  u2]  +  2  a3  vl(0ll)  [ul,  u2]  + 


(1  +  ula)3 

(1+ula)  (a2  v2(0,2>  [ul,  u2]  +2av2(1-0)  [ul,  u2]  +  2  ul  a2  v2a'0)  [ul,  u2]  + 
v2<2'0)  [ul,  u2]  +  2  ul  a  v2|2,0)  [ul,  u2]  +  ul2  a2  v2(2,0)  [ul,  u2] ) )  j} 


■  Take  the  limit  as  v— »0  (or  a — >0)  of  the  NavierStokes  and  divergence  equations  yields 

Limit[nsEq3[ul,  u2] ,  a  ->  0] 

[p(1'0)  [ul,  u2]  , 

p10-1*  [ul,  u2]  +  v2  [ul,  u2]  v2(0'x>  [ul,  u2]  +vl[ul,  u2]  v2 


(1'°>  [ul,  u2]  -  v2(2,0)  [ul,  u2] } 
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‘ Limit [Diverg[{vl[ul,  u2] ,  v2[ul,  u2] }] ,  a  ->  0] 

v2<0'1)  [ul,  u2]  +vl(1’0)  [ul,  u2] 

Thus,  this  implies  the  following  governing  equations  for  the  boundary  layer  flow  in  the  limit  of  small  viscosity: 


Tlr" +  =  0  ( =  vI/2  k(6)  V\2  in  general  case  -  for  this  case*  =  1) 


El4-El  =  n 
d(  m 


where  r  =  t/v112.  These  have  the  same  form  as  Prandtl  limit  equations  for  a  flat  plate.  However,  now  the  boundary 
conditions  are  different.  The  outer  boundary  conditions,  that  is  for  £-><»,  are  found  from  the  matching  conditions  to  the  outer 
layer  solution, which  is  the  inviscjd  solution  around  the  cylinder.  An  analytic  solution  is  known  for  irrotational  flow  which 
provides  the  following  outer  boundary  conditions:  , 


V,  (£,  9)  i-»0,  V2{Z,  9)  —*  -  2  sm9, 

The  inner  boundary  conditions  are  taken  directly  as  the  boundary  conditions  at  the  surface  of  the  cylinder.  Assuming  no-slip 
boundary  conditions  and  control  of  the  normal  velocity  we  have 


Kt(0,  9)  =  c(t,  9),  72(0,0)  =  0, 


pi£,9)  ^  cos20  or  2  cos2  9  or  -2  sin2  9  depending  upon  choice  of  C 
where  c(r,0)  is  the  normal  (control)  velocity  on  the  surface  of  the  cylinder. 


APPENDIX  B:  Source  Listing 


At  the  end  of  this  report  we  have  included  some  source  code  of  the  computer 
implementation  described  in  Section  3.  Printouts  of  the  C++  header  files  of  the  classes 
described  in  the  text  are  listed  and  numbered  individually.  Most  of  these  are  the  headers 
of  the  base  classes  upon  which  the  implementation  is  based.  A  list  of  the  headers  is 
provided  below  along  with  a  brief  description. 

FemLib.h  This  is  the  definition  file  for  the  class  library  FemLib.  It  provides  necessary 
macro  definitions  for  creating  static  and  dynamics  object  libraries.  It  defines  the 
namespace  FemLib  to  which  all  classes  in  FemLib  belong.  It  also  defines  some 
convenience  classes  like  R2  and  the  exception  class  for  the  FemLib  library, 
TSiPemLibError. 

ExtPef.h  Auxiliary  definition  file  containing  extended  type  definitions  (typedefs). 

FeNode2D.h  The  FeNode2D  (actually  named  TSiFeNode2D)  class  definition.  This  class 
is  basically  a  place  marker  for  finite  element  formulations.  It  contains  attributes 
describing  its  various  characteristics  including  its  system  index  and  position  (which  is 
not  really  needed).  Notice  also  that  FeNode2D  maintains  a  list  of  FeElem2D  objects 
to  which  it  is  connected.  The  is  also  an  iterator  class  TSiFeNode2DEIemlterator 
which  is  used  for  iterating  through  all  of  those  objects. 

FeElem2D.h  This  is  the  base  class  for  all  finite  element  types  in  FemLib.  It  is  basically  a 
container  of  FeNode2D  objects  where  the  dependent  variables  are  defined  (or 
sampled  if  you  will).  It  also  includes  an  array  of  vertices,  which  are  R2  objects, 
needed  for  the  differential  geometry  calculations. 

FeQuad2D.h  This  is  an  actual  finite  element  type  in  FemLib  and  is  derived  directly  from 
FeElem2D.  It  implements  all  the  functionality  of  an  arbitrary  quadrilateral  in  the  2D 
plane. 

FeComolex2D.h  The  class  FeComplex2D  is  a  container  of  finite  elements  ( FeElem2D 
objects)  and  nodes  (FeNode2D  objects)  and  represents  a  complex  of  finite  element, 
that  is,  a  geometric  structure  composed  of  polytopes  that  fit  together  "nicely .  It 
provides  methods  for  adding  together  these  objects,  as  well  as  combining  other 
complexes,  to  form  the  complete  object.  The  classes  FeDomain2D  and 
FeBoundary2D  are  derived  from  this  class. 

FeGrid2D.h  Base  class  for  all  finite  element  grids.  A  grid  is  composed  of  an  array  of 
FeDomain2D  objects  and  an  array  of  FeBoundary2D  objects. 

FeBodvInts.h  This  is  a  standalone  class  which  compute  various  integrations  used  in  the 
finite  element  analysis.  The  integrations  are  computed  for  individual  FeElem2D 
objects.  Once  these  integrals  are  computed  for  each  finite  element  they  may  be 
pieced  together  to  form  the  finite  element  system  matrices.  Note  that  a  variety  of 
numeric  quadrature  options  are  possible. 
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Iterators.h  This  is  a  template  file  containing  templates  for  three  different  types  of  iterator 
classes.  One  for  arrays,  one  for  lists,  and  one  for  maps.  It  is  included  because 
iterator  classes  are  frequently  used  in  the  FemLib  library. 
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FemLib.h  3/30/98 


II 

| |  Module  FemLib 

I  I 

\  *========?========= 


*\ 
I  I 
I  I 
I  I 
*/ 


/*  FemLib.h 
* 

*  Macro  Definitions  Module  for  Library  FemLib 

* 

* 

* 

*  Author  :  Christopher  K.  Allen 

*  Copyright  :  August,  1997 

*  Last  Revised  :  September,  1997 

* 


* 


*/ 


/* 

*  Sentry 
*/ 


#ifndef  FemLib_DEFINITION 
#define  FemLib_DEFINITK)N 


*  Create  Library  Importing/Exporting  Macros 
*/ 

#if  defined  (FemLib_EXPORT) 

#  define  FemLiblmpExp  _ declspec  (  dllexport  ) 

#  define  FemLibFunc  _ declspec (  dllexport  ) 

#  define  FemLibClass  _ declspec (  dllexport  ) 

#elif  defined  (FemLib__IMPORT) 

#  define  FemLiblmpExp  _ declspec {  dllimport  ) 

#  define  FemLibClass  _ declspec (  dllimport  ) 

#  if  defined  ( _ FLAT _ ) 

#  define  FemLibFunc  declspec (  dllimport  ) 

#  else 

#  define  FemLibFunc 

#  endif 


#else 

#  define 

#  define 

#  define 


FemLiblmpExp 

FemLibFunc 

FemLibClass 


#endif 


/* 

*  Include  standard  dependencies 
*/ 

// 

//  Avoid  Name  Collisions 

// 

#define  NOMINMAX  //  prevent  min()  and  max()  macros  in  WinDef.h 

#def ine  EXTDEF_NOINDEX  //  prevent  Index  typedef  in  ExtDef.h 


#ifndef  _ AFX_H_ 

#include  <afx.h> 

# endif 

#ifndef  _INC_IO  STREAM 

#include  <iostream.h> 


1 
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#endif 

#ifndef  ExtDef_DEFINITION 
#include  <ExtDef.h> 

#endif  . 

#ifndef  MppExtDef_DEFINITION 
#include  <MppExtDef . h> 

#endif 

#ifndef  MathLib_DEFINITION 
#include  <MathLib\MathLib.h> 
#endif  - 


/*  .  , 

*  Create  the  FemLib  namespace  and  auxiliary  classes 

*/ 


namespace  FemLib  { 

using  namespace  MathLib; 


*  % 


/* 

*  Class  TSiFemLibError 
*/ 


class  FemLibClass  TSiFemLibError  { 
public : 

//  User  Interfaces 
//  Initialization 

TSiFemLibError () ; 

TSiFemLibError (Character*  pErrString) ; 
-TSiFemLibError () ? 

void  SetErrDescr (Character*  pErrString); 

Character*  GetErrDescr ( ) ; 


protected: 

//  Attributes 

TextBuffer  bufDescr;  //  error  description 


} ; 


}? 


#endif 


//  FemLib_DEFINITION 


ExtDef.h  3/30/98 


j  j  Header  ExtDef  | j 


/*  ExtDef.h 

•k 

*  Resource  Module  Containing  Extended  User  Types 

* 

*  Contains  definitions  for  type  extension  to  standard  C++  data 

*  types.  Provides  convenient  precision  modifications  and  indicates 

*  context. 

* 

* 

*  Author  :  Christopher  K.  Allen 

*  Copyright  :  September,  1992 

*  Last  Revised  :  January,  1997 

* 

*/ 


*  Declare  Presence  Identifer 


#ifndef  ExtDef  ^DEFINITION 
#def ine  ExtDef_DEFINITION 


//namespace  ExtDef  { 


//  Alias  basic  data  types 
typedef  char 
typedef  int 
typedef  unsigned 
typedef  double 


Character; 
Integer; 
Cardinal; 
Real ; 


//  New  types 
typedef  void* 
typedef  Character* 
typedef  Character 
typedef  enum  {  False, True  } 


Pointer;  //  generic  pointer 

STRING;  //  beginning  addr  of  string 

TextBuffer [81] ;  //  generic  text  buffer 

Boolean;  //  simple  Boolean,  no  operat's 


//  Situation  specific  types 
typedef  char  Byte; 

typedef  long  int  Filepos; 


//  basic  memory  unit 
//  file  position  type 


#ifndef  EXTDEF__NOINDEX  .  . 

typedef  unsigned  Index;  //  indexing  variable 

#endif 


//}; 


//  Windows  (re) Definitions 
typedef  int  BOOL; 

typedef  unsigned  short  WORD; 
typedef  unsigned  long  DWORD; 

//  Boolean  Variables 

#define  FALSE  0 

#define  TRUE  1 


//  Boolean  variable 
//  machine  word 
//  machine  double  word 


#endif 
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FeNode2D.h  3/30/98 


II 

| |  Definition  TSiFeNode2D 

1 1_ _ _ 


/*  FeNode2D.h 
* 

*  Definition  Module  for  Class  TSiFeNode2D 

* 

*  The  class  TSiFeNode2D  defines  all  the  properties 

*  for  a  grid  point  in  a  finite  element  scheme. 

* 

* 

* 

*  Author  :  Christopher  K.  Allen 

*  Copyright  :  August,  1997  by  Techno  -  Sci  ences ,  Inc. 

*  Last  Modified  :  August,  1997 

* 

*/ 


/* 

*  Declare  Presence  Identifier 


iifndef  FeNode2D_DEFINITI0N 
#def ine  FeNode2D_DEFINITION 


/* 

*  Include  Interface  Dependencies 
*/ 

#ifndef  FemLib_DEFINITION 

^include  <FemLib\FemLib.h> 

#endif 

#ifndef  Iterators_DEFINITION 

#include  <FemLib\Iterators .h> 

#endif 

ttifndef  _ AFXTEMPL_H_ 

iinclude  <afxtempl.h> 

#endif 


/* 

*  Declare  in  namespace  FemLib 
*/ 

namespace  FemLib  { 


/* 

*  Forward  Class  Declaration 
*/ 

class  FemLibClass  TSiFeElem2D; 

class  FemLibClass  TSiFeBody2D; 

class  FemLibClass  TSiFeFace2D; 

class  FemLibClass  TSiFeDomain2D; 

class  FemLibClass  TSiFeBoundary2D; 

class  FemLibClass  TSiFeNode2DElemIterator; 

class  FemLibClass  TSiFeNode2DFaceIterator; 

class  FemLibClass  TSiFeNode2DBodyIterator; 

class  FemLibClass  TSiFeNode2DDomIterator? 

class  FemLibClass  TSiFeNode2DBndIterator; 


1 


==*\ 

II 

II 

II 

==*/ 


FeNode2D.h  3/30/98 


/* 

•*  Class  TSiFeNode2D 
*/' 


class  FemLibClass  TSiFeNode2D  :  public  CObject  { 


DECLARE_SERIAL (  TSiFeNode2D  ) 
public: 

friend  TSiFeNode2DFaceIterator ; 
friend  TSiFeNode2DBody Iterator; 
friend  TSiFeNode2DElemIterator; 
friend  TSiFeNode2DDomIterator; 
friend  TSiFeNode2DBndIterator; 


//  give  face  iterator  access 
//  give  body  iterator  access 
//  give  elem  iterator  access 
//  give  dom  iterator  access 
//  give  bnd  iterator  access 


//  New  Type  Definitions 
enum  GridType  { 
gtGenericPt, 
gtlnternalPt, 
gtBoundaryPt, 
gtCornerPt,  **  * 
nGridTypes, 

1? 

//  enum  NodeType  { 

//  ntGenericNode, 

//  ntCentroidNode, 

//  ntVertexNode, 

//  ntEdgeNode, 

//  nNodeTypes, 

//  > ; 


enum  BndCond  { 
bcNone, 
bcDirichlet, 
bcNeumann, 
bcNatural, 
bcControl, 
bcUserl, 
bcUser2 , 
nBndConds , 

} ; 


// 


typedef  CList<TSiFeElem2D*/  TSiFeElem2D*> 
typedef  CList<TSiFeFace2D*  f  TSiFeFace2D*> 
typedef  CList<TSiFeBody2D* ,  TSiFeBody2D*> 


ElemList; 

FaceList; 

BodyList; 


typedef  CList<TSiFeDomain2D* ,  TSiFeDomain2D*>  DomainList; 

typedef  CList<TSiFeBoundary2D* ,  TSiFeBoundary2D*>  BoundaryList? 


typedef  CArray<RealVector,  RealVector&> 


•  Values Array; 


//  Globals 
static  const  STRING 
static  const  STRING 
static  const  STRING 

static  const  Cardinal 


sGridTypeNames  []  ; 

sNodeTypeNames [] ; 
sBndCondNames  []  ? 

nClassVersion; 


//  string  descriptions 
//  string  descriptions 
//  string  descriptions 

//  pstream  class  version 


public: 

//  User  Functions 
//  Initialization 

TSiFeNode2D()  ;  //  default  constructor 

~TSiFeNode2D() ;  //  destructor 

//  Definition 

void  SetGridType (GridType  gt) ;  //  grid  domain  location 
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//  void  SetNodeType (NodeType  nt) ;  //  set  FE  node  type 

void  SetBndCond (BndCond  be);  //  set  boundary  conditions 

void  SetCoord(R2&  pt) ;  //  set  coordinate 

‘  ,void  SetSys Index (Cardinal  i) ;  //  set  system  node  index 

BOOL  AttachFace (TSiFeFace2D*) ;  //  attach  face  element  to  node 

BOOL  AttachBody (TSiFeBody2D*) ;  //  attach  body  element  to  node 

BOOL  AttachDom(TSiFeDomain2D*) ;  //  attach  domain  to  node 

BOOL  AttachBnd(TSiFeBoundary2D*) ;  //  attach  boundary  to  node 

void  AddValues (RealVectorS  vec) ;  //  add  vector  of  node  values 

//  Assignment 

TSiFeNode2D&  operator= (const  TSiFeNode2D&) ; 

//  Data  Query 

BndCond  GetBndCondO  {  return  bcNode;  }; 

//  NodeType  GetNodeType ( )  {  return  ntNode;  }; 

GridType  GetGridType ( )  {  return  gtNode;  }; 

Cardinal  GetSysIndexO  {  return  iSysIndex;  }; 

R2  GetCoordO  {  return  ptCoord;  }; 

Cardinal  GetNumFaces ( )  {  return  IstFaces .GetCount () ;  }; 

Cardinal  GetNumBod^s ()  {  return  IstBodys .GetCount () ;  }; 

Cardinal  GetNumElems ()  {  return  IstElems .GetCount () ;  };  ✓ 

Cardinal  GetNumDomsO  {  return  IstDoms .GetCount () ;  }; 

Cardinal  GetNumBndsO  {  return  1 stBnds. GetCount () ;  }; 

Cardinal  GetNumValues ()  {  return  arrValues .GetSize () ;  }; 

RealVectorS  GetValues (Cardinal  i) ; 

//  Streams  Support 

virtual  void  StreamDump (ostream&) ;  //  textual  context  dump 

virtual  void  Serialize (CArchiveS  ar) ;  //  pstream  support 


protected: 

//  Local  Data 

GridType  gtNode;  //  node  type  in  domain 

//  NodeType  ntNode;  //  type  of  finite  element  node 

BndCond  bcNode;  //  boundary  conditions  for  node 

Cardinal  iSysIndex;  //  index  into  system 

R2  ptCoord;  //  coordinates  in  R2 

FaceList  IstFaces;  //  list  of  attached  face  elements 

BodyList  IstBodys;  //  list  of  attached  body  elements 

ElemList  IstElems;  //  master  list  of  attached  elements 

BoundaryList  IstBnds;  //  list  of  attached  boundaries 

DomainList  IstDoms;  //  list  of  attached  domains 

ValuesArray  arrValues;  //  value  set  for  node 


protected: 

//  Local  Internal  Functions 
//  Internal  Support 

void  Initialize () ;  //  initialize  data 

BOOL  At tachElem (TSiFeElem2D* ) ;  //  attach  as  general  element 

} ; 


/* 

*  Inline  Functions 
*/ 

// 
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//  Definition 

// 


inline  void  TSiFeNode2D: : SetGridType(GridType  gt) 

,  ,  {  gtNode  =  gt;  }; 

j  M 

//  inline  void  TSiFeNode2D: : SetNodeType (NodeType  nt) 

//  {  ntNode  =  nt;  }; 

inline  void  TSiFeNode2D: : SetBndCond (BndCond  be) 

{  bcNode  =  be;  }; 


inline  void  TSiFeNode2D: : SetSysIndex (Cardinal  i) 
{  iSysIndex  =  i;  }; 

inline  void  TSiFeNode2D: : SetCoord(R2&  pt) 

{  ptCoord  =  pt;  }; 


// 

//  Internal  Support 

// 


// 

//  Stream  Support  *  % 

1 1 


inline  ostream&  operator«  (ostream&  os,  TSiFeNode2D&  e) 
{  e. StreamDump (os) ;  return  os;  }; 


/* 

*  Class  TSiFeNode2DFaceIterator 
*/ 

class  FemLibClass  TSiFeNode2DFaceIterator  : 

public  TSiListIterator<TSiFeNode2D: :FaceList,  TSiFeFace2D> 

{ 

public : 

TSiFeNode2DFaceIterator (TSiFeNode2D*  pNode)  ; 

TSiListIterator<TSiFeNode2D : :FaceList,  TSiFeFace2D> (pNode->lstFaces)  U; 

} ; 


/* 

*  Class  TSiFeNode2DBodyIterator 
*/ 

class  FemLibClass  TSiFeNode2DBodyIterator  : 

public  TSiListIterator<TSiFeNode2D : :BodyList,  TSiFeBody2D> 

{ 

public : 

TSiFeNode2DBodyIterator (TSiFeNode2D*  pNode)  : 

TSiListIterator<TSiFeNode2D : :BodyList,  TSiFeBody2D> (pNode ->lstBodys)  {}; 

}? 


/* 

*  Class  TSiFeNode2DElemIterator 
*/ 


class  FemLibClass  TSiFeNode2DElemIterator  ; 

public  TSiListIterator<TSiFeNode2D : ;ElemList,  TSiFeElem2D> 

{ 

public: 

TSiFeNode2DElemIterator (TSiFeNode2D*  pNode)  : 

TSiLis tIterator<TSiFeNode2D : :ElemListf  TSiFeElem2D> (pNode*>lstElems)  1/ ? 

}? 
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/* 

*  Class  TSiFeNode2DDomIterator 

t/ 

class  FemLibClass  TSiFeNode2DDomIterator  : 

public  TSiListIterator<TSiFeNode2D: :DomainList,  TSiFeDomain2D> 

{ 

public: 

TSiFeNode2DDomIterator (TSiFeNode2D*  pNode)  : 

TSiListIterator<TSiFeNode2D: :DomainList,  TSiFeDomain2D> (pNode ->lstDoms)  {} ; 

} ; 


/* 

*  Class  TSiFeNode2DBndIterator 
*/ 

class  FemLibClass  TSiFeNode2DBndIterator  : 

public  TSiListIterator<TSiFeNode2D: :BoundaryList,  TSiFeBoundary2D> 

{ 

public: 

TSiFeNode2DBndIterator (TSiFeNode2D*  pNode)  : 

TSiListIterator<TSiFeNode2D: :BoundaryList,  TSiFeBoundary2D> (pNode->lstBnds)  {} ; 

};  -  * 


} ;  //  namespace  FemLib 


#endif  //  F eNode2 D_DE F INI T X ON 
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j  j  Definition  TSiFeElem2D  j j 


/*  TSiFeElem2D.h 
* 

*  Definition  Module  for  Class  TSiFeElem2D 

* 

* 

* 

*  Author  :  Christopher  K.  Allen 

*  Copyright  :  August,  1997  by  Techno -Sciences,  Inc. 

*  Last  Revised  :  October,  1997 

* 

*/ 


/* 

*  Include  Sentry 
*/ 

#ifndef  FeElem2D_DEFINITI0N 
idefine  FeElem2D_DEFINITI0N 

-  » 

/* 

*  Include  Interface  Dependencies 
*/ 

#ifndef  FemLib_DEFINITION 
#include  <FemLib\FemLib .h> 

#endif 

#ifndef  FeNode2D_DEFINITION 
#include  <FemLib\FeNode2D.h> 
tendif 

#ifndef  _ AFXTEMP  L_H_ 

#include  <afxtempl.h> 

#endif 


/* 

*  Declare  in  namespace 
*/ 


namespace  FemLib  { 


/* 

*  Forward  Class  Definition 
*/ 

class  FemLibClass  TSiFeElem2DNodeIterator; 
class  FemLibClass  TSiFeElem2DVertlterator; 


/* 

*  Class  TSiFeElem2D 
*/ 


class  FemLibClass  TSiFeElem2D  :  public  CObject  { 
DECLARE_SERIAL (  TSiFeElem2D  ) 


public : 

//  Give  Iterators  Access 

friend  TSiFeElem2DNodeIterator; 
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friend  TSiFeElem2DVertIterator? 


//  New  Types 

enuip  ElemType  {  //  basic  element  types 

enmBasic, 
enmFace, 
enmBody, 
nElemTypes 
}? 

typedef  CArray<R2,  R2&>  FeVertArray; 

typedef  CArray<TSiFeNode2D*,  TSiFeNode2D*>  FeNodeArray; 


//  Global  Data 

static  const  STRING  sElemTypeNames []  ; 

static  const  Cardinal  nClassVersion; 


//  string  descriptors 
//  pstream  class  version 


public : 

//  User  Interfaces 
//  Initialization 

TSiFeElem2D () ; 
virtual  ~TSiFeElem2D ()  ? 

//  Definition  ^  % 

void  SetSysIndex {Cardinal  i Sys Index) ;  //  set  index  into  sys 

virtual  void  SetNode (Cardinal  i,  TSiFeNode2D*  pNode);//  set  node  at  index 
virtual  void  SetVert (Cardinal  iVertex,  R2&  pt) ;  //  set  vertex  coord 


//  Data  Query 
Cardinal  GetNumVerts () 
Cardinal  GetNumNodes ( ) 
Cardinal  GetSysIndex () 


{  return  nVerts;  }; 

{  return  nNodes;  }? 

{  return  i Sys Index;  } ; 


TSiFeNode2D*  GetNode  (Cardinal  i)  ; 
R2&  GetVert (Cardinal  i)  ? 


//  get  ith  node 
//  get  ith  vertex 


//  Event  Hooks 

virtual  void  StreamDump (ostream&  os);  //  textual  content  strmdump 

virtual  void  Serialize (CArchiveS  ar) ;  //  persistent  strms  support 


protected: 


//  Attributes 
ElemType 

enmElemType; 

// 

Cardinal 

nVerts; 

// 

Cardinal 

nNodes ; 

// 

Cardinal 

i Sys Index ; 

// 

FeVertArray 

arrVerts; 

// 

FeNodeArray 

arrNodes ; 

// 

the  basic  element  type 

number  of  element  vertices 
number  of  nodes  on  element 

index  into  FEM  system 

the  element  vertices 
nodes  of  the  element 


protected: 

//  Internal  Support 
void  Initialize ()  ; 

//  For  derived  class 
void  SetElemType  (ElemType  enm)  ; 
void  SetNumNodes  (Cardinal  n)  ; 
void  SetNumVertices  (Cardinal  n)  ; 


//  initialize  parameters 


//  set  basic  element  type 

//  set  number  of  nodes  for  element 

//  set  number  of  vertices  for  elem 


} ; 


/* 

*  Inline  Functions 
*/ 
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// 

//  Definition 

// 

inline  void  TSiFeElem2D: :SetSysIndex (Cardinal  i) 
1  {  iSys Index  =  i;  }; 


// 

//  Internal  Support 

// 

inline  void  TSiFeElem2D: : SetElemType (ElemType  enm) 

{  enmElemType  =  enm;  } ; 


// 

//  Streams  Support 

// 

inline  ostreamS  operator<< (ostreams  os,  TSiFeElem2D&  fe) 
{  fe. StreamDump (os) ;  return  os;  }; 


/*  -  » 

*  Iterator  Classes 


class  FemLibClass  TSiFeElem2DNodeIterator  : 

public  TSiArrayIterator<TSiFeElem2D: :FeNodeArray,  TSiFeNode2D> 

{ 

public : 

TSiFeElem2DNodeIterator (TSiFeElem2D*  pElem)  : 

TSiArrayIterator<TSiFeElem2D: :FeNodeArray,  TSiFeNode2D> (pElem->arrNodes) 

U; 


} ; 


class  FemLibClass  TSiFeElem2DVertIterator  : 

public  TSiArrayIterator<TSiFeElem2D : :  FeVertArray,  R2> 

{ 

public: 

TSiFeElem2DVertIterator (TSiFeElem2D*  pElem)  : 

TSiArrayIterator<TSiFeElem2D: : FeVertArray,  R2> (pElem->arrVerts) 

O; 

} ; 


} ;  //  end  namespace 


#endif  //  FiniteElement_DEFINITION 
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/*===========================*\ 

II  II 

|  j  Definition  TSiFeQuad2D  |  | 


/*  TSiFeQuad2D .h 
* 

*  Definition  Module  for  Class  TSiFeQuad2D 

★ 

* 

* 

*  Author  :  Christopher  K.  Allen 

*  Copyright  :  August,  1997  by  Techno -Sciences,  Inc. 

*  Last  Revised  :  August,  1997 

* 

*/ 


/* 

*  Include  Sentry 
*/ 

#ifndef  FeQuad2D_DEFINITI0N 
#define  FeQuad2D_DEFINITION 

-  » 

/* 

*  Include  Interface  Dependencies 
*/ 

#ifndef  FeBody2D_DEFINITION 
#include  <FemLib\FeBody2D.h> 
#endif 


/* 

*  Declare  in  namespace 
*/ 


namespace  FemLib  { 

//  using  namespace  ExtDef; 

//  using  namespace  MppExtDef; 


/* 

*  Forward  Class  Definition 
*/ 

class  FemLibClass  TSiFeFace2D; 


/* 

*  Class  TSiFeQuad2D 
*/ 

class  FemLibClass  TSiFeQuad2D  s  public  TSiFeBody2D  { 

DECLARE_SERIAL (  TSiFeQuad2D  ) 

public: 

//  New  Types 

typedef  TSiFeBody2D: : Interpolation  Interpolation; 

//  Global  Data 

static  const  Cardinal  nClassVersion;  //  pstream  class  version 


public: 

//  User  Interfaces 
//  Initialization 
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TSiFeQuad2D () ; 

~TSiFeQuad2D ()  ; 

virtual  TSiFeFace2D*  MakeFace (Cardinal  i) ;  //  make  ith  body  face 

.  virtual  void  MakeAllFaces () ;  //  make  all  body  faces 

//  FEM  Requirements 

virtual  void  Setlnterpolation (Interpolation  intp) ; 

virtual  Real  ShapeFunc (Cardinal  iNode,  R2&  ptLoc) ;  //  function  value 

virtual  R2  ShapeGrad (Cardinal  iNode,  R2&  ptLoc);  //  function  deriv 

virtual  R2  GetCoordValue (R2&  ptLoc);  //  get  coord  from  local  coord 

virtual  R2x2  GetCoordDeriv(R2&  ptLoc);  //  get  coord  derivative 

virtual  Real  GetCoord Jacob (R2&  ptLoc);  //  get  coord  xform  Jacobian 

virtual  BOOL  Membership (R2&  ptGbl) ;  //  pt  in  element? 

//  Event  Hooks 

virtual  void  StreamDump  (os  t  reams  os); 
virtual  void  Serialize (CArchive&  ar)  ; 


protected; 

//  Attributes 


protected: 

//  Internal  Support 
void  Initialize ()  ; 

//  Shape  functions  for  interpolation  scheme 

Real  ShapeFuncConst (Cardinal  iNode,  Real  r,  Real  s) ; 

Real  ShapeFuncBiLin (Cardinal  iNode,  Real  r.  Real  s) ; 

Real  ShapeFunc Seren (Cardinal  iNode,  Real  r,  Real  s) ; 

Real  ShapeFuncQuad (Cardinal  iNode,  Real  r,  Real  s) ; 

R2  ShapeGradConst (Cardinal  iNode,  Real  r,  Real  s)  ; 

R2  ShapeGradBiLin (Cardinal  iNode,  Real  r,  Real  s) ; 

R2  ShapeGradSeren (Cardinal  iNode,  Real  r,  Real  s) ; 

R2  ShapeGradQuad (Cardinal  iNode,  Real  r.  Real  s)  ; 

//  Used  by  Membership ( ) 

BOOL  InsideSegment (R2&  ptGbl,  R2&  ptSegVl,  R2&  ptSegV2) ; 

}  ? 


/* 

*  Inline  Functions 
*/ 


// 

//  FEM  Requirements 

// 


// 

//  Coordinate  Functions 

// 


// 

//  Streams  Support 

// 

inline  ostream&  operator« (ostream&  os,  TSiFeQuad2D&  fe) 
{  fe. StreamDump (os) ;  return  os;  }; 


} ;  //  end  namespace 
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#endif  //  FeQuad2D_DEFINITION 


-  * 


FeComplex2D.h  3/30/98 


I  1  II 

|  |  Definition  TSiFeComplex2D  |  | 


/*  FeComplex2D.h 

k 

*  Definition  Module  for  Class  TSiFeComplex2D 

* 

* 

*  TSiFeComplex2D  is  the  container  class  for  a  finite  element  system 

*  representation  in  2  dimensions.  It  contains  TSiFeNode2D  objects 

*  and  TSiFeElem2D  Objects  which  define  the  system. 

* 

*  NOTE : 

*  A  TSiFeComplex2D  object  will  not  allow  duplicate  nodes  or  elements. 

*  The  members  AddNodeO  and  AddElemO  will  return  FALSE  if  such  an 

*  attempt  is  made. 


k 

k 

k 

k 

k 

k 

*/ 


Author 
Copyright 
Last  Revised 


:  Christopher  K.  Allen 

:  September,  1997  by  Techno- Sciences ,  Inc. 
:  October,  1997 


-  i 


/* 

*  Include  Sentry 
*/ 

#ifndef  FeComplex2D_DEFINITION 
#define  FeComplex2D_DEFINITION 


/*  ,  .  y 

*  Include  Interface  Dependencies 
*/ 

#ifndef  FeElem2D__DEFINITI0N 
tinclude  <FemLib\FeElem2D ,h> 

#endif 

#ifndef  FeNode2D__DEFINITION 
#include  <FemLib\FeNode2D.h> 

#endif 

#ifndef  Iterators_DEFINITION 
#include  <FemLib\Iterators .h> 

#endif 

#ifndef  _ AFXTEMPL__H _ 

#include  <afxtempl.h> 

#endif 


/* 

*  Declare  in  namespace 
*/ 


namespace  FemLib  { 


/* 

*  Forward  class  definitions 
*/ 

class  FemLibClass  TSiFeComplex2DNodeIterator; 
class  FemLibClass  TSiFeComplex2DElemIterator; 
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/* 

*  Class  TSiFeComplex2D 
*/ 

class  FemLibClass  TSiFeComplex2D  :  public  CObject  { 
DECLARE_SERIAL(  TSiFeComplex2D  ) 
public: 

//  Give  the  iterators  access 
friend  TSiFeComplex2DNodeIterator; 

friend  TSiFeComplex2DElemIterator; 


//  New  Types 

typedef  CArray<TSiFeNode2D* ,  TSiFeNode2D*>  FeNodeArray; 
typedef  CArray<TSiFeElem2D* ,  TSiFeElem2D*>  FeElemArray; 


//  Global  Data 

static  const  Cardinal  nClassVersion;  //  pstream  class  version 


public : 

//  User  Interfaces 
//  Initialization 

TSiFeComplex2D ( ) ; 
virtual  ~TSiFeComplex2D ( ) ; 


//  Definition 

void  SetHomogeneous (BOOL  bHomog) ; 

BOOL  AddNode (TSiFeNode2D*  pNode) ; 
BOOL  AddElem(TSiFeElem2D*  pElem) ; 


//  all  elements  same  size 

//  add  node  to  grid 
//  add  element  to  grid 


BOOL  AddComplex(TSiFeComplex2D*  p)  ;  //  add  objects  of  complex 


void  RemoveAl 1 { )  ; 
void  DestroyAll ()  ? 


//  remove  all  objects 
//  destroy  entire  grid 


//  Data  Query 

BOOL  GetHomogeneous ( ) 

Cardinal  GetNumElems ( ) 

Cardinal  GetNumNodes ( ) 

TSiFeNode2D*  GetNode (Cardinal  i) ; 
TSiFeElem2D*  GetElem (Cardinal  i) ; 


{  return  bHomog ;  } ; 

{  return  nElems;  }; 
{  return  nNodes ;  } ; 


Index  GetNodeSys  Index  ( )  ? 


//  Event  Hooks 

virtual  void  StreamDump (ostream&  os); 
virtual  void  Serialize (CArchive&  ar) ? 


protected: 

//  Attributes 

BOOL  bHomog;  //  is  homogeneous  grid  ? 


Cardinal  nElems; 

Cardinal  nNodes; 


//  number  of  element  vertices 
//  number  of  nodes  on  element 


FeElemArray 

FeNodeArray 


arrElems; 
arrNodes  ? 


//  array  of  elements 
//  array  of  nodes 


protected: 

//  Internal  Support 
void  Ini tialize ( ) ; 

BOOL  FindNode (TSiFeNode2D*  p)  ;  //  true  if  node  is  in  complex 

BOOL  FindElem (TSiFeElem2D*  p) ;  //  true  if  element  is  in  complex 


2 


FeComplex2D.h  3/30/98 


/* 

*  Inline  Functions 

*/ 

// 

//  Definition 

// 

inline  void  TSiFeComplex2D: : SetHomogeneous (BOOL  b) 
{  bHomog  =  b ;  } ; 


// 

//  Streams  Support 

// 

inline  ostream&  opera tor«  (os tream&  os,  TSiFeComplex2D&  fe) 
{  fe.StreamDump(os)  ;  return  os;  }; 


/* 

*  Class  TSiFeComplex2DNodeIterator 
*/ 

class  FemLibClass  TSiFeComplex2DNodeIterator  : 

public  TSiArrayIterator<TSiFeComplex2D: :FeNodeArray,  TSiFeNode2D> 

{ 

public: 

TSiFeComplex2DNodeIterator (TSiFeComplex2D*  pGrid)  : 

TSiArrayIterator<TSiFeComplex2D: :  FeNodeArray,  TSiFeNode2D>  (pGrid ->arrNodes)  {} ; 


} ; 


/* 

*  Class  TSiFeComplex2DElemIterator 
*/ 

class  FemLibClass  TSiFeComplex2DElemIterator  : 

public  TSiArrayIterator<TSiFeComplex2D: : FeElemArray,  TSiFeElem2D> 

{ 

public: 

TSiFeComplex2DElemIterator (TSiFeComplex2D*  pGrid)  :  ,  ,  .  r, 

TSiArrayIterator<TSiFeComplex2D:  : FeElemArray,  TSiFeElem2D>  (pGnd->arrElems)  {}  ; 


} ; 


} ;  //  end  namespace 

#endif  //  NsSystem2D_DEFINITI0N 
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M  II 

j  j  Definition  TSiFeGrid2D  |  | 


/*  TSiFeGrid2D.h 
* 

*  Definition  Module  for  Class  TSiFeGrid2D 

* 

* 

*  TSiFeGrid2D  is  the  container  class  for  a  finite  element  system 

*  representation  in  2  dimensions.  It  contains  TSiFeNode2D  objects 

*  and  TSiFeElem2D  Objects  which  define  the  system. 

* 

*  NOTE : 

*  A  TSiFeGrid2D  object  OWNS  all  the  objects  it  contains.  Therefore, 

*  once  it  goes  out  of  scope  all  its  objects  are  destroyed. 

* 

* 

* 

*  Author  :  Christopher  K.  Allen 

*  Copyright  :  September,  1997  by  Techno* Sciences ,  Inc. 

*  Last  Revised  :  September,  1997 

* 

*/ 

-  % 

/* 

*  Include  Sentry 
*/ 


#ifndef  FeGrid2D_DEFINITI0N 
#def ine  FeGrid2D_DEFINITI0N 


*  Include  Interface  Dependencies 
*/ 

#ifndef  FeBoundary2D_ DEFINITION 
#include  <FemLib\FeBoundary2D . h> 
#endif 

tfifndef  FeOb j ect2D_DEFINITION 
#include  <FemLib\FeObj ect2D.h> 
#endif 

#ifndef  FeDomain2D_DEFINITI0N 
#include  <FemLib\FeDomain2D . h> 
#endif 

iifndef  IteratorsJDEFINITION 
#include  <FemLib\ Iterators .h> 
#endif 

#ifndef  _ AFXTEMPL_H_ 

#include  <afxtempl .h> 

#endif 


/* 

*  Declare  in  namespace 
*/ 


namespace  FemLib  { 


/* 

*  Forward  class  definitions 
*/ 
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class  FemLibClass  TSiFeGrid2DNodeIterator; 
class  FemLibClass  TSiFeGrid2DElemIterator; 


/* 

*  Class  TSiFeGrid2D 
*/ 

class  FemLibClass  TSiFeGrid2D  :  public  CObject  { 

DECLARERS  ERIAL  (  TSiFeGrid2D  ) 
public: 

//  Give  the  iterators  access 
friend  TSiFeGrid2DNodeIterator ? 

friend  TSiFeGrid2DElemIterator ; 

//  New  Types 

typedef  CArray<TSiFeDomain2D* ,  TSiFeDomain2D*>  DomArray; 
typedef  CArray<TSiFeBoundary2D*,  TSiFeBoundary2D*>  BndArray; 
typedef  CArray<TSiFeOb j  ect2D* ,  TSiFeOb j  ect2D*>  Ob  j  Array; 


//  Global  Data 

static  const  Cardinal  nClassVersion?  //  pstream  class  version 

-  % 


public : 

//  User  Interfaces 
//  Initialization 

TSiFeGrid2D ( ) ; 
virtual  ~TSiFeGrid2D ( ) ; 


//  Definition 

void  AddDomain  (TSiFeDomain2D*  pDom) ;  //  add  a  domain  to  grid 

void  AddBoundary  (TSiFeBoundary2D*  pBnd)  ;  //  add  domain  boundary  to  grid 

void  AddObject (TSiFeOb ject2D*  pObj);  //  add  embedded  object  to  grid 


void  RemoveDomain (Cardinal  iDom) ; 
void  RemoveBoundary (Cardinal  iBnd); 
void  RemoveObject (Cardinal  iOb j ) ? 


//  remove  domain  at  index 
//  remove  boundary  at  index 
//  remove  boundary  at  index 


void  DestroyAll () ? 


//  destroy  entire  grid 


void  TagCornerNodes  ()  ; 


//  identify  corner  nodes 


return  cpxMaster .GetNumNodes () ;  }; 
return  cpxMaster .GetNumElems () ;  }; 

return  arrDoms .GetSize () ;  }; 
return  arrBnds . GetSize ( ) ?  } ; 
return  arrObjs. GetSize () ;  }; 

TSiFeBoundary2D*  GetBoundary (Cardinal  iBnd); 

TSiFeOb j ect2D*  GetObj ect (Cardinal  iBnd) ; 

TSiFeDomain2D*  GetDomain (Cardinal  iDom  =  0)  ; 

TSiFeComplex2D&  GetMasterComplex () ; 

//  Event  Hooks 

virtual  void  StreamDump (ostream&  os); 
virtual  void  Serialize (CArchive&  ar) ; 


//  Data  Query 

Cardinal  GetNumNodes ( )  { 

Cardinal  GetNumElems ( )  { 

Cardinal  GetNumDomains ( )  { 

Cardinal  GetNumBoundaries ()  { 

Cardinal  GetNumObjects ()  { 


protected: 


//  Attributes 
TSiFeComplex2D 

cpxMaster; 

// 

the  master  complex 

DomArray 

arrDoms ; 

// 

array  of 

domains 

BndArray 

arrBnds ; 

// 

array  of 

domain  boundaries 

Ob j Array 

arrObjs; 

// 

array  of 

free  boundaries 

protected: 
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//  Internal  Support 
void  InitializeO ; 
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} ; 


/* 

*  Inline  Functions 

*/ 

// 

//  Definition 

// 


// 

//  Streams  Support 

// 

inline  ostreams  operator« (ostreamS  os,  TSiFeGrid2D&  fe) 
{  fe. StreamDump (os) ;  return  os;  }; 


/* 

*  Class  TSiFeGrid2dNodpIterator 
*/ 

class  FemLibClass  TSiFeGrid2DNodeIterator  : 
public  TSiFeComplex2DNodeIterator 

{ 

public: 

TSiFeGrid2DNodeIterator (TSiFeGrid2D*  pGrid)  : 

TSiFeComplex2DNodeIterator (  & (pGrid- >cpxMaster)  )  {}  ; 

}  ; 


/* 

*  Class  TSiFeGrid2DElemIterator 
*/ 


class  FemLibClass  TSiFeGrid2DElemIterator  : 
public  TSiFeComplex2DElemIterator 

{ 

public: 

TSiFeGrid2DElemIterator (TSiFeGrid2D*  pGrid)  : 

TSiFeComplex2DElemIterator (  & (pGrid ->cpxMas ter)  )  {}  ; 

} ; 


} ;  //  end  namespace 


iendif  //  FeGrid2D_DEFINITION 


3 


FeBodyInts2D.h  3/30/98 


ii  '  ~  ii 

|  j  Definition  TSiFeBodyXnts2D  |  1 


/*  FeBodyInts2D .h 
* 

*  Definition  Module  for  Class  TSiFeBodyInts2D 

* 

* 

*  Computes  the  shape  function  integrals  in  the  finite  element 

*  solution  of  the  viscous,  incompressible  Navier- Stokes 

*  equations . 

* 

* 

*  Author  :  Christopher  K.  Allen 

*  Copyright  :  August,  1997  by  Techno- Sciences ,  Inc. 

*  Last  Revised  :  October,  1997 

* 

*/ 


/* 

*  Include  Sentry 
*/ 

-  f 

#ifndef  FeBodyInts2D_DEFINITION 
#def ine  FeBodyInts2D_DEFINITION 


/* 

*  Include  Interface  Dependencies 
*/ 

#ifndef  FeBody2D_DEFINITION 
#include  <FemLib\FeBody2D.h> 
#endif 

#ifndef  _ AFXTEMPL.H _ 

#include  <afxtempl ,h> 

#endif 


/* 

*  Declare  in  namespace 
*/ 


namespace  FemLib  { 

//  using  namespace  ExtDef; 

//  using  namespace  MppExtDef; 


/* 

*  Class  TSiFeBodyInts2D 
*/ 


class  FemLibClass  TSiFeBodyInts2D  :  public  CObject  { 
DECLARE_SERIAL {  TSiFeBodyInts2D  ) 


public: 

//  New  Types 
enum  Quadrature  { 
enmNone , 
enmSimpson, 
enmBode, 
enmGaussian3 , 
enmGaussian4 , 
enmGaussian5 , 
enmExt Simpson, 


//  quadrature  not  set 
//  exact  for  3rd  order  poly 
//  exact  for  5th  order  poly 
//  good  for  4th  order  poly 
//  good  for  6th  order  poly 
//  good  for  8th  order  poly 
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nQuadratures 
}  ? 

typedef  CArray<R2,  R2&>  QuadPtArray; 
typedef  CArray<Real,  Real>  QuadWtArray; 


public : 

//  Global  Data 

static  const  STRING  sQuadratureNames [] ;  //  quadrature  descr.'s 

static  const  Real  f SimpsonPts  [] ?  //  bode  quad  points 

static  const  Real  f SimpsonWts  [] ;  //  bode  quad  weights 

static  const  Real  fBodePts [] ;  //  bode  quad  points 

static  const  Real  fBodeWts [] ;  //  bode  quad  weights 

static  const  Real  fGaussianPts3  [] ;  //  gauss3  quad  points 

static  const  Real  fGaussianWts3 [] ?  //  gauss3  quad  weights 

static  const  Real  fGaussianPts4 [] ;  //  gauss4  quad  points 

static  const  Real  fGaussianWts4 [] ;  //  gauss4  quad  weights 

static  const  Real  fGaussianPts5  [] ;  //  gaussS  quad  points 

static  const  Real  fGaussianWts5 [] ;  //  gauss5  quad  weights 

static  const  Cardinal  nClassVersion;  //  pstream  class  version 


public : 

//  User  Interfaces^  % 

//  Initialization 

TSiFeBodyInts2D ( )  ; 

~TSiFeBodyInts2D ()  ; 

//  FEM  Configuration 

void  SetQuadrature (Quadrature  qd) ;  //  set  num  quadrature 

//  FEM  Computations 

BOOL  CompBasisMatrix (TSiFeBody2D*  pFe) ;  //  comp  stiffness  matrix 

BOOL  CompGradMatrix (TSiFeBody2D*  pFe) ;  //  comp  viscous  mat 

BOOL  CompAdvecTensors (TSiFeBody2D*  pFe) ;  / /  comp  x  advection  ten 

BOOL  CompDivMatrices  (TSiFeBody2D*  pFeVel,  //  comp  x  divergence  mat 

TSiFeBody2D*  pFePres) ; 

//  Data  Query 

RealMatrixS  GetBasisMatrixO  {  return  matBasis;  }; 

RealMatrix&  GetGradMatrix ( )  {  return  matGrad;  }; 

RealMatrix&  GetDivXMatrix ()  {  return  matDivX;  }; 

RealMatrix&  GetDivYMatrix ()  {  return  matDivY;  }; 

RealTensor&  GetAdvecXTensor ()  {  return  tenAdvecX;  }; 

RealTensor &  GetAdvecYTensor ()  {  return  tenAdvecY;  }; 

//  Event  Hooks 

virtual  void  StreamDump (ostream&  os); 

virtual  void  Serialize (CArchiveS  ar) ? 


protected: 

//  Attributes 
//  Configuration 

Quadrature  enmQuadrature;  //  numeric  quadrature  type 

//  The  integral  matrices 

RealMatrix  matBasis; 

RealMatrix  matGrad; 

RealMatrix  matDivX; 

RealMatrix  matDivY; 

RealTensor  tenAdvecX; 

RealTensor  tenAdvecY; 

//  Quadrature  variables 

Cardinal  nQuadOrder;  //  polynomial  order  of  quadrature 

Cardinal  nQuadPts;  //  number  of  quadrature  pts  total 

Real*  pQuadPts ;  //  array  of  standard  quadrature  pts 

Real*  pQuadWts;  //  array  of  standard  quadrature  wts 

QuadPtArray  arrQuadPts;  //  quadrature  points 

QuadWtArray  arrQuadWts;  //  quadrature  weights 


//  velocity  basis  matrix 
//  viscosity  matrix 
//  divergence  matrix  in  x  dir 
//  divergence  matrix  in  y  dir 
//  advection  tensor  for  x  vel 
//  advection  tensor  for  y  vel 
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protected: 


//  Internal  Support 
void  Initialize ()  ; 

1 

}; 

i 

void  Ini tQuadrature ( )  ; 

//  init  numeric  quad 

/* 

★ 

Inline  Functions 

*/ 

// 

// 

// 

FEM  Requirements 

// 

// 

// 

Streams  Support 

inline  ostream&  operator«  (ostream&  os, 

TSiFeBodyInts2D&  fe) 

{ 

f e. StreamDump (os) >  ^return  os?  }; 

} ; 

//  end  namespace 

#endif 

//  FiniteElementJDEFINITION 
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/*=  \ 
j  j  Definition  TSilterators  |  j 


/*  iterators. h 
* 

*  Definition  Module  for  Class  TSilterators 

* 

* 

* 

*  Author  :  Christopher  K.  Allen 

*  Copyright  :  January,  1997  by  Techno -Sciences,  Inc. 

*  Last  Revised  :  April,  1997 

* 

*/ 


/* 

*  Include  Sentry 
*/ 


#ifndef  Iterators_DEFINITION 
#define  Iterators_DEFINITION 


*  Include  Interface  Dependencies 
*/ 


#ifndef  _ AFX_H _ 

#include  <afx.h> 
#endif 

#ifndef  _ AFXTEMPL_H _ 

tinclude  <afxtempl .h> 
#endif 


/* 

*  Declare  in  FemLib  Namespace 
*/ 

//namespace  FemLib  { 

//  using  namespace  ExtDef? 


/ 


*  * 
* 

* 

* 

*/ 


Class  TSiArraylterator 


template  <class  Array,  class  Object> 
class  TSiArraylterator  { 

public: 

//  User  Interfaces 
//  Initialization 

TSiArraylterator (Array&  rArray); 
^TSiArraylterator  () ; 

void  Restart (); 

//  Data  Query 
Object*  GetCurrent  0 

Object*  operator++(int)  ; 

operator  BOOLO  ? 


//  restart  the  iterator 

{  return  pCurrObject;  }; 

//  return  current  Object  and  advance 
//  still  valid  Objects's  ? 


protected: 
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//  Attributes 

Array& 

rArray; 

Cardinal 

pos; 

Object* 

pCurrObject; 

}■; 

i 

// 

// 

// 

Member  Functions 

//  the  attached  map  class 
//  current  position  in  map 
//  the  current  argument 


template  <class  Array,  class  Object> 

TSiArrayl terator<Array ,  Object>: :  TSiArrayl terator(Array&  __rArray)  : 
rArray  (__rArray) 

{ 

pos  =  0; 

pCurrObject  =  NULL; 

} 


template  <class  Array,  class  Object> 

TSiArrayl  ter  at  or  <Ar  ray ,  Ob  j  ect> : :  -TSiArrayl  t  era  tor  ( ) 


{} 


template  <class  Array,  class  Object> 
void  TSiArrayl terator<A*ray,  Object>:  : Restart  () 
{ 

pos  =  0 ; 

pCurrObject  =  NULL; 

} 


template  <class  Array,  class  Object> 

Object*  TSiArrayl terator<Array,  Object>: :operator++ (int) 

if  (  pos  <  (Cardinal) rArray. GetSize ()  ) 

pCurrObject  =  ( Ob j ect*) rArray .GetAt (  pos++  ); 
else 

pCurrObject  =  NULL; 
return  pCurrObject; 

} 


template  <class  Array,  class  Object> 

TSiArrayl terator<Array,  Ob j ect> :: operator  BOOL() 

return  (  pos  <  rArray. GetSize ( )  ); 

} 


* 


*  Class  TSiListlterator 

* 

*/ 


***★ 


template  <class  List,  class  Object> 
class  TSiListlterator  { 

public: 

//  User  Interfaces 
//  Initialization 
TSiListlterator (Lists  _rList) ; 

-TSiListlterator () ; 

void  Restart ();  //  restart  the  iterator 

//  Data  Query 

Object*  GetCurrentO  {  return  pCurrObject;  }; 

Object*  operator** (int) ?  //  return  current  Assoc  and  advance 
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operator  BOOL ( ) ; 


//  still  valid  Assoc* s  ? 


protected: 


. 

//  Attributes 

Lists 

rList; 

POSITION 

pos; 

}; 

Object* 

pCurrObj ect; 

// 

// 

// 

Member' Functions 

//  the  attached  map  class 
//  current  position  in  map 
//  the  current  argument 


template  <class  List,  class  Object> 

TSiListIterator<List,  Object>: :TSiListIterator (Lists  _rList)  : 
rList  (_rList) 

{ 

pos  =  rList .GetHeadPosition  ()  ; 
pCurrObject  =  NULL; 

} ; 


template  <class  List,  class  Object> 
TSiListIterator<List,  Gbje ct>: :~TSiListIterator () 

O; 


template  <class  List,  class  Object> 
void  TSiListIterator<List,  Object>: :Restart () 
{ 

pos  =  rList .GetHeadPosition () ; 
pCurrObject  =  NULL; 

} ; 


template  <class  List,  class  Object> 

Object*  TSiListIterator<List,  Object>: :operator++ (int) 
{ 

if  (pos) 

pCurrObject  =  (Ob j ect*) rList .GetNext (  pos  ); 
else 

pCurrObject  =  NULL; 
return  pCurrObject; 

}? 


template  <class  List,  class  Object> 
TSiListIterator<List,  Ob j ect> :: operator  BOOLO 
{ 

return  (pos  !=  NULL); 

}; 


/**********★★******************************************* 

★ 

*  Class  TSiMapIterator 

★ 

*/ 

template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 
class  TSiMapIterator  { 

public : 

//  User  Interfaces 
//  Initialization 

TSiMapIterator (Map&  rMap); 

-TSiMapIterator ()  ; 

void  Restart ()?  //  restart  the  iterator 
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//  Data  Query 

Assoc*  GetCurrentO  {  return  pCurrAssoc;  }; 


Assoc*  operator** (int) ; 
i  operator  BOOLO; 


//  return  current  Assoc  and  advance 
//  still  valid  Assoc 's  ? 


protected: 


//  Attributes 
Maps 

rMap  ; 

// 

the 

attached  map  class 

POSITION 

pos; 

//  current  position  m  map 

} ; 

Assoc* 

pCurrAssoc; 

// 

the 

current  argument 

// 

// 

Member  Functions 

// 

template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 
TSiMapIterator<Map,  Key,  Assoc,  Map_Assoc> :  :TSiMapIterator  (Map&  _rMap)  : 
rMap  (_rMap) 

{ 

pos  =  rMap.GetStartPositionO; 
pCurrAssoc  =  NULL; 

} ;  -  * 


template  <class  Map,  class  Key, 
TS iMapI t era  tor <Map ,  Key,  Assoc, 
{}; 


class  Assoc,  class  Map_Assoc> 
Map_Assoc>: : -TSiMapIterator  () 


template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 
void  TSiMapIterator<Map,  Key,  Assoc,  Map_Assoc> ::  Res  tart  () 


}; 


pos  =  rMap.GetStartPositionO; 

pCurrAssoc  =  NULL; 


template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 

Assoc*  TS iMapI terator<Map ,  Key,  Assoc,  Map_Assoc>: :operator++ (int) 

{ 

Key  temp; 

Map_Assoc*  pMapAssoc; 

if  (pos) 

{ 

rMap. GetNext Assoc (  pos,  temp,  pMapAssoc  ); 
pCurrAssoc  =  (Assoc*) pMapAssoc; 

} 

else 

pCurrAssoc  =  NULL; 
return  pCurrAssoc; 

} ; 


template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 
TSiMapIterator<Map,  Key,  Assoc,  Map_Assoc>: : operator  BOOLO 
{ 

return  (pos  !=  NULL) ; 

}; 


* 

*  Class  TSiMapKeylterator 

* 
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*/ 


template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 
class  TSiMapKeylterator  { 


public: 

//  User  Interfaces 
//  Initialization 

TSiMapKeylterator (Map&  rMap) ? 
-TSiMapKeylterator () ? 

void  Restart (); 

//  Data  Query 
Key  GetCurrent ( ) 

Key  operator++ (int)  ? 

operator  B00L0  ; 


//  restart  the  iterator 
{  return  key;  } ? 

//  return  current  Assoc  and  advance 
//  still  valid  Assoc' s  ? 


protected: 


//  Attributes 
Map& 

rMap  ; 

// 

the  attached  map  class 

POSITION 

pos? 

// 

current  position  in  map 

}  /• 

Key 

key; 

*  % 

// 

the  current  key 

// 

// 

Member  Functions 

// 

template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 
TSiMapKeyIterator<Map,  Key,  Assoc,  Map_Assoc>:  : TSiMapKeylterator  (Map&  _ rMap) 
rMap  LrMap) 

pos  =  rMap . GetStartPosition { ) ; 

//  key  =  NULL; 

}? 


template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 
TSiMapKeyIterator<Map,  Key,  Assoc,  Map_Assoc> :: -TSiMapKeylterator () 


{}? 


template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 
void  TSiMapKeyIterator<Map,  Key,  Assoc,  Map__Assoc>:  :  Restart  () 

pos  =  rMap .GetStartPosition () ? 

//  key  =  NULL; 


} ; 


template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc>  ^ 
Key  TS iMapKeyl t erator<Map ,  Key,  Assoc,  Map_Assoc>:  :operator++(mt) 
{ 

Map_As  s  oc  *  pMap As  s  oc ; 
if  (pos) 

rMap . GetNext Assoc (  pos,  key,  pMapAssoc  ); 

} 

//  else 

//  pCurr Assoc  =  NULL; 

return  key; 

}? 


template  <class  Map,  class  Key,  class  Assoc,  class  Map_Assoc> 
TS  iMapKeyl  ter  a  tor<Map,  Key,  Assoc,  Map_Assoc>: :  operator  B00L() 
{ 

return  (pos  1=  NULL) ; 

} ; 
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//}? 

#endif 


//  namespace  FemLib 
//  I terator s_DEFINITION 
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1  Introduction 


In  this  short  note  we  present  a  possible  technique  for  designing  feedback  control  systems 
for  viscous  fluid  flow.  The  technique  is  based  on  the  analysis  of  spatial  modes  generated 
by  objects  in  the  fluid  domain  and  uses  the  group  properties  of  the  trigonometric 
expansion  functions.  We  start  with  the  two-dimensional  vorticity  transport  equations  and 
expand  the  field  variables  in  terms  of  trigonometric  exponentials.  The  time  dependent 
equations  for  the  expansion  coefficients  are  then  derived  from  the  governing  equations. 
The  structure  of  these  equations  leads  to  a  strategy  for  controlling  the  fluid  motion. 


2  The  Vorticity  Transport  Equations 


From  the  two-dimensional,  incompressible  Navier-Stokes  equations  we  may  derive  the 
vorticity  transport  equations.  These  equations  describe  the  time  evolution  and  spatial 
distribution  of  the  fluid’s  vorticity  co.  This  behavior  is  coupled  to  that  of  the  fluid’s  stream 
function  vy  as  shown  below. 

d>-vV2co  +  J(co,ii/)  =  0  ^ 

-V2!//  =co 


2.1 


where  the  over  dot  indicates  time  differentiation,  v  is  the  kinematic  viscosity  of  the  fluid, 
and  J(oa,\\i)  is  the  Jacobian  determinant  of  co  and  y  (i.e.,  J((i\\\i)=d(d\\\j)/d(x,y) ).  The  scalar 
functions  co  andy  are  related  to  the 
vector  valued  velocity  field  v  by  the 
following:  / 

Ly - 

<o  =az  •  V x v,  v=Vx(azyr).  (2) 


Therefore,  we  can  recover  the  fluid’s 
velocity  field  from  the  above  once 
system  (1 )  is  solved. 

The  Basis  Functions 


4  x 


Consider  the  region  shown  in  Figure 
1  depicting  a  rectangular  fluid  domain 

in  cartesian  coordinates  (x,y).  The  •  Figure  1 :  Fluid  domain 

dimensions  are  Lx  and  Ly, 


respectively.  We  let  ©  and  \\i  be  elements  of  the  function  space  spanned  by  the  set  of 
basis  functions 


(3) 


where 


n  =  (/t|,K,)eZxZ 
x  =  (xi,x!)  =  (xj)€R! 


2  n  Ik 
nx — ,n2  — 
1  Lr  2  Lv 


(4) 


Remarks: 

1)  The  basis  functions  e‘k"1  are  the  eigenfunctions  of  the  Laplacian  operator  V2  with 
eigenvalues  —  |kn[2  =  -(jfc2  +  k'l ).  This  makes  inversion  of  the  Laplacian  operator 
particularly  easy. 

2)  The  basis  ^,k"  *}  also  forms  a  group  under  multiplication  which  is  isomorphic  to  the 
group  Z  ©  Z  with  the  isomorphism 

e'k"  *  h-»(«„n2).  (5) 

3)  The  function  space  that  {?ik"  *}  generates  is  then  isomorphic  to  the  group  ring 
c(z  ©  z) .  These  properties  are  used  when  computing  the  Jacobian  J(co,y).  The 
nonlinearities  generated  by  the  Jacobian  term  are  represented  by  multiplication  within 
the  ring. 


2J2  Computation  of  the  Expansion  Coefficients 

The  vorticity  ©  and  stream  function  vy  may  be  expanded  in  terms  of  the  above  basis 
functions. 


-ko  ±2, 

a>(x,y;t)  =  ^©n(r)e'k”x  and  y{x,y\t)  =  ■  *  (6) 

n=(«1,rt2)=-«  n=(«j  ,«2 

where 

ffln(0  =  -r=r  Jj(o(x,y;t),e-lkt'tdxdy, 

4LxLy  [0,LxH0,Ly]  ~ 

y.(0=  ,  * —  jj\i/(x,y;t),e~‘K  ldxdy. 

[0.ij]x[0,lj,] 

The  spatial  variations  are  represented  by  the  basis  functions  while  the  time  variation  is 
contained  in  the  expansion  coefficients  an(t)  and  yn(f). 
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We  insert  these  expansions  into  the  vorticity  transport  equations  to  get  first-order 
differential  equations  for  the  expansion  coefficients.  Note  that  from  the  second  of  Eqs.  (1) 
we  may  express  the  in  terms  of  the  by  the  orthogonality  of  the  basis  functions 

*  We  have 


Using  this  fact  we  need  only  consider  the  vorticity  coefficients.  Substituting  the 
expansions  (6)  into  (1)  and  using  relation  (8)  yields  (after  carrying  out  the  group 
multiplications) 


CO 


,  ,  k,  k,  -  k„  k,, 

<  +v(k\  +kni)(0„  -  X  —  2  ,2--J-<P,<»t=0 

s+t=n  Ktx  +  Kt2 


(9) 


We  have  one  such 
equation’for  every  mode 
coefficient  co„.  (Recall 
that  n,  s,  and  t  are  all  in 
the  space  Z  ©  Z .)  The 
second  term  on  the  LHS 
of  (9)  represents  the 
self-damping  caused  by 
the  fluid  viscosity  while 
the  third  term  (the 
summation)  constitutes 
the  dynamics  of  the 
coupling  between 
modes  in  the  fluid. 

Referring  to  Figure  2  we 
see  schematically  how 
low  order  modes  can 
excite  higher  order  modes  through  this  process. 


Figure  2:  Schematic  of  mode  coupling 


3  Controlling  Fluid  Flow 


Now  consider  an  object  embedded  in 
the  fluid  domain,  say  an  airfoil.  The 
situation  is  depicted  in  Figure  3.  Also, 
assume  that  there  is  some  sort  of 
sensor/actuator  assembly  located  on 
the  airfoil  that  we  are  able  to  control. 
We  wish  to  determine  which  spatial 
modes  the  airfoil  can  excite  and, 
conversely,  which  modes  the 
actuators  can  excite.  Knowing  this 
behavior  we  hope  to  design  a  control 
law  which  suppresses  modes  causing 
the  most  difficulty. 


•  Figure  3:  Airfoil  in  fluid  domain 
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3.1  Determining  Mode  Excitation 


Assume  that  the  boundary  of  the  airfoil,  y,  can  be  parameterized  according  to  the 
following: 


Y : [0,1] [0, Lx]x [0, Ly ] 

n-^(x(r),y(r)) 


(10) 


Now  the  set  of  restrictions  je*"'*  |y|  forms  a  basis  for  the  space  of  complex  analytic 

functions  on  the  boundary  y,  however,  this  basis  is  not  orthogonal.  We  can  use  this 
condition  to  compute  which  modes  the  airfoil  boundary  will  intrinsically  excite.  Let  the 
function  coo  defined  on  y  represent  the  no-slip  boundary  conditions  on  the  airfoil 
(determination  of  coo  is  not  straight  forward  but  assume  for  the  sake  of  argument  that  it 
exists — boundary  layer  theory  will  probably  be  needed).  Then  we  have  the  requirement 

.  i  a>(x,  y;  O  =  o0  on  y.  (11) 

This  requirement  determines  which  modes  are  intrinsically  excited  by  the  airfoil.  We  may 
compute  this  excitement  (approximately)  by  multiplying  both  sides  of  Eq.  (11)  by  e~,k“'x 
then  integrating  over  the  boundary  y.  Truncating  at  some  /V=(A/i,/V2)  yields  the  matrix 
vector  equation 


\e*",eik'x)i  - 

•*,<?*’  *}  •••  (eik"  x,e‘k*  ^ 

1'  'r  ) 

where 

1 

//.  >  /.•>!  (13) 


and 

(«*■■*, ©o)  =jffl0(r)e-v"(x(r,-v<r,W.  (14) 

0 

Equation  (12)  can  be  inverted  using  standard  techniques  to  determine  the  excited  mode 
strengths. 

Now  let  the  function  aa(()  defined  on  y  represent  the  action  of  the  actuators  on  the  airfoil. 
The  boundary  conditions  for  the  actuators  require  that 

a(x, y;t)  =  a, 2(t)  on  y  (15) 

for  all  time  t.  With  a  similar  argument  as  before,  the  modes  and  their  strengths  that  the 
actuators  excite  at  a  particular  time  f  is  computed  by  the  following  equation: 
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(e'k'\e'k'^  ... 

•••  (e,ltN  V*N  l) 


'  <S)r  ^ 


jf®»v 


Now  that  we  know  the  modes  that  we  have  control  over  and  their  relative  strengths,  along 
with  the  modes  that  are  excited  naturally,  we  should  be  able  to  design  a  controller  using 
Eq.  (9)  as  a  guide  for  the  mode  dynamics. 

3.2  Other  Application 

It  is  possible  to  apply  these  same  techniques  to  the  Navier-Stokes  equations  directly  in 
their  primitive  variable  form.  The  results  in  that  case  will  not,  however,  be  as  compact  as 
they  are  here.  Yet  they  may  be  more  conducive  to  practical  application.  This  may  be 
especially  true  if  we  consider  the  restricted  case  of  the  time-varying  boundary  layer 
equations  developed  in  the  Phase  I  final  report. 

-  % 
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